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ABSTRACT 


In  the  Chemical,  Biological,  Radiological -Nuclear  and  Explosives  Research  and  Technology  Initiative 
(CRTI)  Project  09-509TD,  a  persistent,  highly  realistic,  game-based  synthetic  environment  (SE)  for 
incident  commanders  and  first  responders  to  conduct  individual  and/or  team  training  and  collaboration 
(at  the  tactical  and/or  strategic  levels)  to  address  deliberate  or  accidental  Chemical,  Biological, 
Radiological,  Nuclear  and  Explosives  (CBRNE)  releases  in  an  urban  environment  has  been  proposed. 
This  will  be  achieved  by  integrating  a  state-of-the-science  physics-based  urban  flow/dispersion 
modeling  system  and  a  high-fidelity  building-aware  Geographic  Information  System  (GIS) 
representation  of  a  real  cityscape  (The  City  of  Calgary)  with  a  proven  game-based  simulation  engine 
developed  by  3DIntemet  Inc.  WATCFD  is  responsible  for  the  provision  of  hazardous  plume  entities  to 
the  game-based  simulation  engine.  The  objective  of  this  report  is  to  provide  a  technical  description  of 
modeling  and  simulation  of  the  highly  disturbed  building-induced  flow  field  in  downtown  Calgary  and 
the  transport  and  dispersion  of  hazardous  agents  released  into  this  complex  flow  field  using  the 
computational  fluid  dynamics  (CFD)  model  urbanSTREAM,  which  consists  of 

1 .  generation  of  a  structured  computational  grid  from  a  building  database  for  an  urban  landscape  in 
the  form  of  an  Arc  View  Shape  file  for  the  City  of  Calgary; 

2.  generation  of  highly  disturbed  building-aware  flow  fields  in  the  urban  environment  for  various 
inflow  conditions; 

3.  development  of  a  heavy-gas  modeling  capability  and  implementing  it  in  the  dispersion  module 
urbanEU; 

4.  simulation  of  dispersion  for  each  of  the  CBRNE  hazard  scenarios  described  in  the  portfolio  and 
archiving  the  concentration  field  as  a  function  of  space  and  time; 

5.  conversion  of  the  concentration  field  datasets  from  various  simulations  into  a  “compact  data” 
format  compatible  for  utilization  in  the  game-based  simulation  environment  developed  by 
3DIntemet  Inc.; 

6.  provision  of  input/support  for  presentations  in  the  2012  annual  Public  Security  S&T  Summer 
Symposium 

It  is  envisaged  that  the  development  of  the  game-based  virtual  reality  CBRNE  training  environment 
developed  in  the  present  project  will  provide  the  highest  quality  training  to  incident  commanders  and 
first  responders,  allowing  the  end  user  to  experience  and  respond  to  realistic  high-risk  scenarios 
involving  CBRNE  hazards  in  a  complex  cityscape  in  a  completely  safe  and  controlled  environment. 
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EXECUTIVE  SUMMARY 


First  Responder  Immersive  Training  Simulation  Environment 

Fue-Sang  Lien1’2 3 4 5,  Kun-Jung  Hsieh1  and  Hua  Ji1 

Waterloo  CFD  Engineering  Consulting  Inc.  (WATCFD) 

2:  University  of  Waterloo 

Background:  Considerable  efforts  have  been  made  to  strengthen  Canada’s  preparedness  for,  prevention 
of,  response  to,  and  recovery  from  Chemical,  Biological,  Radiological,  Nuclear  and  Explosives 
(CBRNE)  threats  to  public  safety  and  security.  An  important  part  of  this  effort  involves  CBRNE  training 
for  first  responders  in  order  to  provide  them  with  the  knowledge,  competence  and  confidence  to  deal 
with  a  wide  spectrum  of  CBRNE  hazards.  Presently,  scenario-based  training  in  Canada  is  limited  to  first 
responders’  hands-on  training  (including  live-agent  training)  at  Defence  R&D  Canada,  Suffield 
Research  Centre  and  the  Canadian  Emergency  Management  College  (CEMC).  However,  this  training 
does  not  utilize  simulated  urban  environments  with  realistic  hazard  scenarios.  Consequently,  there  is  a 
significant  capability  gap  in  providing  realistic  CBRNE  training  for  first  responders.  In  this  project,  a 
game-based  synthetic  environment  will  be  developed  to  address  this  capability  gap.  3DIntemet  Inc.  will 
integrate  their  unique  simulation  engine  with  a  downwind  hazard  dispersion  model,  a  Geographic 
Information  System  (GIS)  representation  of  Calgary,  and  expert  first  responder  input  to  produce  the 
overall  simulation  environment  for  training.  The  development  of  a  portfolio  of  CBRNE  hazard  scenarios 
(e.g.,  a  derailment  and  explosion  involving  volatile  chlorine  tanker  in  downtown  Calgary,  etc.)  will  be 
led  by  LuomaTech  Inc.,  and  all  partners  in  the  project  will  be  involved  in  the  portfolio  development. 
WATCFD  will  utilize  building-resolving  urban  flow/dispersion  models  to  provide  hazardous  plume 
entities  for  the  game-based  simulation  engine  to  be  developed  by  3DIntemet  Inc.  This  training 
environment  will  enable  first  responders  to  train  together  more  frequently  against  a  wider  spectrum  of 
realistic  CBRNE  hazard  scenarios. 

Results:  This  report  provides  a  technical  description  of  the  modeling  and  simulation  of  the  highly 
disturbed  building-induced  flow  field  in  downtown  Calgary  and  the  transport  and  dispersion  of 
hazardous  agents  released  into  this  complex  flow  field  using  the  computational  fluid  dynamics  model 
urbanSTREAM  developed  by  WATCFD  which  consists  of 

1 .  generation  of  a  structured  computational  grid  from  a  building  database  for  an  urban  landscape  in 
the  form  of  an  Arc  View  Shape  file  for  the  City  of  Calgary; 

2.  generation  of  the  highly  disturbed  building-aware  flow  fields  in  the  urban  environment  for 
various  inflow  conditions; 

3.  development  of  a  heavy-gas  modeling  capability  and  implementing  it  in  the  dispersion  module 
urbanEU; 

4.  simulation  of  dispersion  for  each  of  the  CBRNE  hazard  scenarios  described  in  the  portfolio  and 
archiving  the  concentration  field  as  a  function  of  space  and  time; 

5.  conversion  of  the  concentration  field  data  sets  from  the  various  simulations  into  a  “compact 
data”  format  compatible  for  utilization  in  the  game -based  simulation  environment  developed  by 
3DIntemet  Inc.; 


iii 


6.  provision  of  input/support  for  presentations  in  the  2012  annual  Public  Security  S&T  Summer 
Symposium 

Significance:  It  is  envisaged  that  the  development  of  the  game-based  virtual  reality  CBRNE  training 
environment  developed  in  the  present  project  will  provide  the  highest  quality  training  to  incident 
commanders  and  first  responders,  allowing  the  end  user  to  experience  and  respond  to  realistic  high-risk 
scenarios  involving  CBRNE  hazards  in  a  complex  cityscape  in  a  completely  safe  and  controlled 
environment. 
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1.  INTRODUCTION 


Since  the  terrorist  attacks  of  September  11,  2001,  and  the  subsequent  establishment  of  the 
Chemical,  Biological,  Radiological-Nuclear,  and  Explosives  (CBRNE)  Research  and 
Technology  Initiative  (CRTI),  considerable  efforts  have  been  made  to  strengthen  Canada's 
preparedness  for,  prevention  of,  response  to,  and  recovery  from  CBRNE  threats  to  public  safety 
and  security.  An  important  part  of  this  effort  involves  CBRNE  training  for  first  responders  in 
order  to  provide  them  with  the  knowledge,  competence  and  confidence  to  deal  with  a  wide 
spectrum  of  CBRNE  incidents.  At  the  present  time,  scenario-based  training  is  limited  to  first 
responders’  hands-on  training  at  Defence  R&D  Canada,  Suffield  Research  Centre  (Counter 
Terrorism  Technology  Centre  -  CTTC)  and  other  courses  offered  by  the  Canadian  Emergency 
Management  College  (CEMC).  A  new  CRTI  project  (09-509TD)  has  been  funded  to  augment 
our  current  capabilities  to  train  first  responders  to  respond  to  CBRNE  incidents. 

The  objective  of  this  project  is  to  develop  a  highly  realistic,  game -based  simulation  environment 
for  first  responders  to  conduct  individual  and/or  team  training  and  collaborations  in  addressing 
deliberate  or  accidental  releases  of  CBRNE  materials  in  a  real  cityscape.  This  will  be  achieved 
by  integrating  a  physics-based  urban  flow  and  dispersion  modeling  system  developed  in  previous 
CRTI  projects  and  a  high-fidelity  GIS  representation  of  an  urban  environment  with  a  proven 
game-based  simulation  engine  to  produce  a  simulation  environment  for  CBRNE  training / 
education  that  is  realistic  and  engaging. 

Under  the  current  project,  the  following  five  tasks  have  been  completed  by  WATCFD,  which 
will  be  addressed  in  detail  in  the  present  report. 

1.  Generate  an  appropriate  structured  computational  grid  from  a  building  database  in  the 
form  of  an  Arc  View  Shapefile  for  the  City  of  Calgary  for  the  use  of  the  subsequent 
simulation  of  the  flow  and  dispersion  of  materials  released  into  this  urban  environment; 

2.  For  a  given  set  of  inflow  conditions  (mean  wind  speed,  wind  direction,  and  turbulence 
intensities),  simulate  the  highly  disturbed  building-aware  flow  field  in  the  urban 
environment  and  generate  building-aware  flow  fields  for  various  inflow  conditions; 

3.  For  a  given  portfolio  of  CBRNE  hazard  scenarios,  and  using  the  disturbed  flow  fields 
generated  in  Task  2  above,  simulate  the  dispersion  for  each  of  these  scenarios  in  the 
urban  environment  and  archive  the  concentration  field  as  a  function  of  space  and  time; 

4.  Convert  the  concentration  field  datasets  from  the  various  simulations  in  a  format 
compatible  for  utilization  in  the  game-based  simulation  environment  developed  by 
3DIntemet  Inc.  (http://www.3dintemet.com/); 

5.  Provide  input/support  for  presentations  and/or  technical  papers  prepared  for  various 
scientific  meetings  and  symposia  (including  the  annual  Public  Security  S&T  Summer 
Symposium). 

In  addition,  a  heavy-gas  (or  dense-gas)  modeling  capability  has  been  also  developed  in  this 
project,  which  is  relevant  to  Scenario  2  (to  be  discussed  in  Sections  3  and  4),  in  which  the 
derailment  of  a  rail  car  containing  chlorine  is  considered.  Chlorine  is  a  heavier-than-air  (dense) 
gas. 
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2.  TECHNICAL  DESCRIPTION 


2.1  Overview  of  the  Integrated  Urban  Flow  and  Dispersion  Modeling  System 

The  urban  modeling  system  developed  from  a  previous  CRTI  Project  02-0093RD  includes  five 
main  modules:  urbanGRID,  urbanSTREAM,  urbanEU,  urbanAEU,  and  urbanPOST.  These 
modules  and  how  they  interface  with  each  other  and  with  other  project  components  have  been 
summarized  in  Figure  2  of  the  report  by  Yee  et  al.  (2007).  As  stated  in  this  report,  urbanGRID 
imports  building  information  encoded  in  Environmental  Systems  Research  Institute  (ESRI) 
Shapefiles  and  uses  this  data  to  generate  a  structured  grid  over  a  user-selected  computational 
domain  in  a  given  cityscape.  Furthermore,  urbanGRID  imports  three-dimensional  meteorological 
fields  (e.g.,  mean  wind,  turbulence  kinetic  energy,  etc.)  provided  by  urban  GEM  LAM1  and  uses 
this  information  to  provide  inflow  boundary  conditions  for  the  urban  microscale  flow  model: 
urbanSTREAM.  The  structured  grid  and  inflow  boundary  conditions  provided  by  urbanGRID  are 
used  as  inputs  by  urbanSTREAM  to  simulate  the  flows  around  and  within  the  complex 
geometries  of  buildings  in  the  cityscape,  and  to  provide  the  high-resolution  wind  and  turbulence 
fields  for  dispersion  modeling.  The  two  Eulerian  dispersion  models  urbanEU  (source-oriented) 
and  urbanAEU  (receptor-oriented)  use  the  wind  and  turbulence  fields  provided  by 
urbanSTREAM  to  simulate  the  dispersion  of  contaminants  in  the  urban  domain.  Finally, 
urbanPOST  is  the  post-processing  module  to  process  the  primary  output  files  from 
urbanSTREAM,  which  can  then  be  visualized  by  third  party  visualization  programs  such  as 
Tecplot2. 

To  consider  effect  from  complex  terrain  on  a  flow  as  required  in  the  present  project, 
urbanSTREAM  applies  curvilinear  coordinates,  allowing  for  the  use  of  a  body-fitted  (or 
boundary-fitted)  mesh  to  include  the  topographic  capability.  The  governing  equations  written  in 
curvilinear  coordinates  can  be  converted  from  those  written  in  Cartesian  coordinates  based  on  a 
one-to-one  (locally  invertible)  mapping  (or  transformation).  For  example,  all  the  governing 
partial  differential  equations  under  Cartesian  coordinates  can  be  written  in  a  generic  form  as 
follows, 


d(f>  Surf 

dt  dxj 


d0_ 

dxn 


(1) 


for  a  scalar  </> ,  which  can  be  identified  with  any  of  the  mean  Cartesian  velocity  components  u,  , 

turbulence  kinetic  energy  k,  viscous  dissipation  s ,  mean  concentration  C,  or  influence  function 
C*.  Transforming  the  above  equation  from  3-D  Cartesian  coordinates  to  3-D  curvilinear 
coordinates,  represented  by  ^'(y  =  1,2,3  and  ^  =rj ,  =  C,  ),  produces  a  partial 

different  equation  in  the  following  generic  form, 


1  urban  GEM  LAM  is  a  prognostic  mesoscale  model  with  an  urban  parameterization  developed  by 
Environment  Canada  for  CRTI  Project  02-0093RD. 

2  Tecplot  (http://www.tecplot.com/)  is  a  commercially-available  CFD  visualization  software. 
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where  m,n  =  1,2,3 ,  are  indices  following  Einstein  summation  convention,  J  is  the  Jacobian 
matrix  of  the  transformation,  that  is, 
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and  UJ  ( j  =  1, 2, 3 )  are  the  contravariant  velocities,  expressed  as  follows, 


U2=V  =  Uur/  +  vr/  +  wr/)  =  j(u^-  +  v^-  +  w^- 

v  x  y  z}  \  dx  dy  dz ; 
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in  which  u,v,w  denote  the  mean  Cartesian  velocity  components. 
The  fourth-order  tensor  qt is  defined  as 

€i=p]j: 


(4) 

(5) 

(6) 


(V) 


where  /?/  (i,  j  =  1,2,3)  are  the  elements  of  the  following  inverse  Jacobian  matrix, 
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(8) 


The  transformation  of  the  source  term  Sy/>  is  very  tedious  and  its  end  result  is  different  for  each 
governing  equation.  Due  to  space  constraints,  detailed  expressions  of  transformed  are  not 
included  in  this  section,  and  interested  readers  can  refer  to  Lien  (1992). 
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2.2  Inflow  Boundary  Condition 

Instead  of  taking  inflow  velocity  profiles  from  GEM  LAM  solutions,  the  velocity  profiles  at  the 
inlet  (or,  inflow  boundary)  are  approximated  by  the  following  power  law: 


u(z) 


u(z) 


z  y 

—  for  z  <  400  m 

vioj 


(400f 

10 


for  z  >  400  m 


(9) 


based  on  the  assumption  that  p  «  0.3 ,  and  uw  ~  5.556  m/s  (or  ul0  ~  20  km/h  ),  where  «l0  is  the 
wind  speed  at  10  m  above  the  ground.  The  exponent  g«0.3  in  equation  (9)  is  a  good 
approximation  for  wind  speeds  in  suburbs.  The  profile  of  turbulence  kinetic  energy  at  the  inflow 
boundary  is  expressed  by  using  the  following  two-layer  model: 
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(11) 


In  equations  (10)  and  (11),  hc  «  40  m  is  the  average  building  height  in  the  inner  region  shown  in 
Figures  7  and  8  later,  and  hBL  « 1000  m  is  the  atmospheric  boundary-layer  thickness.  The 
maximum  turbulence  kinetic  energy  kmax  »  0.4  m2  /  s2  is  estimated  from  the  IOP9  test  problem 
(Allwine  et  al.,  2004)  in  the  Joint  Urban  2003  (JU2003)  experiment  (https://ju2003- 
dpg.dpg.army.mil/').  and  kmm  «  O.Olu2,  at  z  »  hBL  .  Finally,  the  turbulence  dissipation  rate,  e,  is 
approximated  by 

,3/2 

^  =  ^71 - - -  (12) 

Cfl3  4  min(x*z,  hBL  /  3) 

where  k  =  0.42  is  the  von  Karman  constant,  and  Cu  =0.09  for  the  standard  k  —  s  turbulence 
closure  model. 

The  resulting  system  of  partial  differential  equations  was  solved  numerically  using  a  collocated, 
finite-volume  method  and  implemented  in  the  urbanSTREAM  code  (Yee  et  al.,  2007;  Lien  et  al., 
2010).  Diffusive  volume-face  fluxes  were  discretized  using  a  second-order  accurate  central 
differencing  scheme  (CDS).  Advective  volume-face  fluxes  were  approximated  using  a  second- 
order  accurate  Upstream  Monotonic  Interpolation  for  Scalar  Transport  (UMIST)  scheme  (Lien 
and  Leschziner,  1994).  The  transient  term  was  discretized  using  a  fully  implicit,  second-order 
accurate  three-time-level  method  described  in  Ferziger  and  Peric  (2002).  The  Semi-Implicit 
Method  for  Pressure-Linked  Equations  (SIMPLE)  algorithm  (Patankar  and  Spalding,  1972)  was 
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used  to  determine  the  (gauge)  pressure.  A  nonlinear  interpolation  scheme  (Rhie  and  Chow,  1983) 
was  used  to  interpolate  the  cell  face  velocities  from  the  adjacent  nodal  velocities  at  the  cell 
centers  in  order  to  prevent  checkerboard  oscillations  from  developing  in  the  pressure  field. 


2.3  Inclusion  of  Heavy  Gas  Capability 

2.3.1  Mathematical  formulation  and  numerical  framework 

The  dispersion  of  a  heavier-than-air  (dense)  gas  cloud  is  predicted  by  solving  the  three- 
dimensional  conservation  equations  for  mean  (ensemble-averaged)  quantities  in  a  turbulent  flow 
field.  This  involves  solving  the  conservation  equations  of  mass,  momentum,  and  species 
concentration  for  an  isothermal  (constant  temperature)  fluid  flow  based  on  the  Reynolds- 
averaged  Navier-Stokes  (RANS)  approach: 


0£+A(pjj)=o, 

8t  dXjK  j) 

8t V  dxj  '  7  dxj  dxj 


f  8u^ 

\  uxj  y 


8  / _ \  8 

ato)+a^»A)-a^ 


pd S 

.  dxn 
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dx 


-  —  {puV^  +  ip- P»)g,, 
-[pufk],  k  =  l,2,...,N, 


(13) 

(14) 

(15) 


where  ui  is  the  ensemble-mean  (or,  Reynolds-averaged)  velocity  in  the  x-direction  with  i  =  1,2, 
3  representing  the  x,  y,  and  z  directions,  respectively;  t  is  the  time;  p  is  the  mean  (gauge) 
pressure;  p  is  the  dynamic  molecular  viscosity  of  the  fluid;  ck  is  the  mass  fraction 
(concentration)  of  the  A>th  scalar  (or  species)  in  the  mixture;  D  is  the  molecular  diffusivity  of  the 
scalar;  g/  is  the  gravitational  acceleration  in  the  xrdirection;  p  is  the  density  of  the  fluid;  and,  p() 
is  a  reference  density  (chosen  herein  as  the  density  of  air  at  sea  level  with  a  representative  value 
of  p{)  =  1 .25  kg  m~3).  Summation  is  implied  by  repeated  indices.  Reynolds-averaging  of  the 
momentum  transport  equation  (14)  and  the  advection-diffusion  equation  (15)  gives  rise, 
respectively,  to  the  Reynolds  stresses  which  are  defined  as  the  tensor  up.  and  the  turbulent 

scalar  fluxes  which  are  defined  as  the  vector  UjCk  (k  =  1,2,. .  ,,N).  These  equations,  along  with  an 
equation  of  state  that  uses  the  ideal  gas  law  approximation  for  the  density,  namely 


P  = 


P 


N 


T^Rkck 


P 

RS’ 


(16) 


are  the  main  governing  equations.  Here,  P  is  the  absolute  pressure,  T  is  the  absolute  temperature 
(assumed  to  be  a  constant  for  the  isothermal  atmosphere  considered  herein),  and  Rk  =  G/ Mk  is 
the  specific  gas  constant  for  species  k  (G  =  8.314472xl03  J  K'1  kg-mof1  is  the  universal  gas 
constant  and  Mk  is  the  molecular  weight  of  species  k).  Also,  in  equation  (16),  Rt  is  the  specific 
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gas  constant  for  the  mixture  consisting  of  N  species.  It  is  noted  that  the  gauge  pressure  p  in 
equation  (14)  is  the  pressure  deviation  from  an  adiabatic  atmosphere  in  hydrostatic  equilibrium 
with  corresponding  density  p0 ,  with  the  result  that  the  absolute  pressure  in  equation  (16)  is 

related  to  the  gauge  pressure  in  equation  (14)  as  P  =  p  +  P0  ,  where  P0  is  some  nominal 
atmospheric  pressure  level  in  the  adiabatic  atmosphere  (e.g.,  at  the  ground  surface). 

If  the  density  variations  are  not  too  large,  we  can  apply  the  Boussinesq  approximation  to 
equations  (13),  (14)  and  (15).  The  main  step  in  developing  the  Boussinesq  approximation  is  to 
treat  the  density  as  a  constant  (which  can  be  taken  as  the  reference  density  p() )  in  the  unsteady, 
convective  and  diffusive  terms,  while  retaining  the  density  variations  only  in  the  buoyancy 
(gravity)  term  (p- p{) )gi  in  the  mean  momentum  transport  equation.  This  approach  was  used  by 

Perdikaris  and  Mayinger  (1993)  and  Scargiali  et  al.  (2005)  to  simulate  dense  gas  dispersion  over 
a  topographically  complex  terrain.  However,  the  approximation  for  the  buoyancy  term  given  in 
these  two  papers  is  not  correctly  stated.  In  consequence,  we  provide  a  simple  derivation  of  the 
Boussinesq  approximation  of  (p  -  p(}  )gj  in  the  context  of  the  dispersion  of  a  heavier-than-air  gas 
released  into  the  atmosphere.  Towards  this  purpose,  consider  a  Taylor  series  expansion  of  the 
density  p  =  p(P,T,cl,...,cN)  to  first  order  about  a  reference  state  (denoted  using  a  subscript  0) 
given  by 


P  =  Po  + 


dp 

dP 
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P0)+  — 

'  dT 


(T~T»)+t§ 

I  k= 1  UCk 


(Ci  Ck, 0  )  ’ 


(17) 


with  the  partial  derivatives  in  equation  (17)  evaluated  at  the  reference  state  ( P0,T0,c10,...,cNfi ) 

and  p0  =  p(P0,T0,clfi,...,cNQ) .  For  an  isothermal  atmosphere,  the  third  term  on  the  right-hand 

side  of  equation  (17)  vanishes  identically.  Inserting  the  ideal  gas  law  of  equation  (16)  into 
equation  (17),  the  density  perturbations  (from  the  reference  state)  can  be  expressed  as  follows: 


P~Po  

Po  Po 


(P-Pq) 

Po 


~Ck,o) 
k= 1 
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Po 
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(18) 


where  Sp  =  (p-p0)  and  SP  =  (P-P0)  =  p  are  used  to  denote  the  density  and  pressure 
perturbations,  respectively;  and, 


1  dp 

Po  Pck  0 


Pk 


i= 1 


(19) 


is  the  volumetric  expansion  coefficient  arising  from  perturbations  in  the  concentration  of  species 
k. 
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To  simplify  the  approximation  in  equation  (18),  we  consider  order-of-magnitude  estimates  for 
SP/P0  in  comparison  to  Sp/ p{) .  To  proceed,  it  is  assumed  that  the  pressure  gradient  term  in  the 

mean  momentum  transport  equation  (14)  is  of  the  same  order  of  magnitude  as  the  buoyancy  term 
and,  with  the  buoyancy  force  acting  in  the  y  (i  =  2)  direction,  this  implies  within  the  Boussinesq 
approximation  that 


1  dp 

1  dSP 

t 

i 

Po  £y 

po  £y 

5 

Po 

g 

Po 

(20) 


An  estimate  of  the  magnitude  of  the  term  on  the  left-hand- side  of  equation  (20)  [obtained  by 
using  the  ideal  gas  law  of  equation  (16)  for  the  reference  state]  gives 


I  dSP  ~  RJ ;  |£P| 

Po  £,y  po  h 


(21) 


where  ly  is  a  representative  or  characteristic  spatial  scale  of  the  circulations  in  the  atmospheric 
surface  layer  (and  typically,  ly  =  0(100)  m  in  an  adiabatic  atmosphere).  Insertion  of  the  order-of- 
magnitude  estimate  of  equation  (21)  into  equation  (20)  leads  to 

\8P\  _  lyg  \Sp\  _  /,g  \Sp\  _  lygp0  \Sp\ ' 

P0  R  J{)  p0  P0/p0  p0  P0  p0  ‘ 


To  proceed  further,  we  need  an  estimate  for  Pq  which  can  be  obtained  as  follows.  The  reference 
pressure  Po  at  the  Earth’s  surface  is  determined  for  an  atmosphere  in  hydrostatic  equilibrium,  so 
Po  ~  pogHo  where  H0  is  the  effective  height  in  the  atmosphere  where  the  pressure  is  equal  to  zero 
(approximately).  Ho  is  estimated  to  be  about  Ho  =  0(8000)  m  (Pielke,  2002).  Now,  using  this 
estimate  for  P0  in  equation  (22),  it  follows  that 


M  lySPo  M  _  ly8Po  M  _  h  \SP\ 
P0  P0  Po  Po&H o  Po  Ho  Po 


(23) 


from  which  it  is  seen  that  \SP\j P0  « \5p\j p{)  because  ly/H0  «1.  From  this  observation,  it  is 
now  evident  that  the  density  perturbations  in  equation  (18)  are  well  approximated  as 


P-Po  =  Sp 
Po  Po 


(24) 


Finally,  for  the  special  case  of  a  two-component  mixture  consisting  of  air  (a,  or  species  k  =  1) 
and  a  contaminant  species  (g,  or  species  k  =  2)  [chlorine,  Freon- 12  or  some  other  contaminant  in 
the  cases  considered  herein],  the  density  perturbation  of  equation  (24)  reduces  to 
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(25) 


P~Po 

Po 


=  ac , 


on  noting  ca  =  1  —  c  =  1  —  c,  a  =  (p  - pa)/ pg  ,  and  assuming  that  the  concentrations  of  air  and 
contaminant  gas  in  the  reference  state  are  ca  0  =  1  and  c  0  =  1  -  ca  0  =  0  .  Using  the  density 
perturbation  given  by  equation  (25),  the  buoyancy  term  can  be  approximated  as  follows: 

(p~Po)gi=Po<Kgi.  (26) 

For  the  simulation  of  the  dense  gas  dispersion  scenarios  in  this  report,  equation  (26)  will  be  used 
to  model  the  buoyancy  term  in  the  mean  momentum  equation  within  the  framework  of  the 
Boussinesq  approximation.  It  is  noted  that  this  approximation  for  the  buoyancy  term  differs  from 
that  used  in  Scargiali  et  al.  (2005)  in  that  the  quantity  a  involves  normalization  by  pa  rather  than 
pg  (as  used  in  the  current  formulation).  Finally,  for  a  two-component  mixture  consisting  of  air 
and  a  contaminant  gas,  we  need  only  consider  a  transport  equation  for  the  contaminant  species 
concentration  c  because  the  concentration  of  the  air  in  the  mixture  is  determined  from  (1  -  c) . 


Within  the  framework  of  the  standard  high-Reynolds-number  k-s  model,  the  Reynolds  stresses 
up  t  and  the  turbulent  scalar  fluxes  up  ,  required  to  close  the  transport  equations  for  the  mean 

momentum  and  mean  concentration  of  the  contaminant  species,  are  modeled  using  the 
Boussinesq  eddy-viscosity  approximation  (which  should  not  be  confused  with  the  Boussinesq 
approximation  for  buoyant-driven  flows)  and  the  simple  gradient  diffusion  hypothesis, 
respectively,  as  follows: 


—  2  ( dut  du  ^ 

up ,  =  —  ko,.  -  v.  — -  H - - 

1  ]  3  9  \dxj  dxtJ 

tt  v  dc 

UJC  = - 

ac  dxj 


(27) 

(28) 


where  vt  =  C/t  k2  je  is  the  kinematic  (turbulent)  eddy  viscosity,  Cu  =  0.09  is  a  model  constant,  crc 

=  0.63  is  the  turbulent  Schmidt  number  for  the  scalar,  k  is  the  turbulence  kinetic  energy,  s  is  the 
dissipation  rate  of  k,  and  Sy  is  the  Kronecker  delta  function.  The  modeled  transport  equations  for 
k  and  s  in  the  standard  k-s  model  with  buoyancy  terms  assume  the  following  form  within  the 
context  of  the  Boussinesq  approximation: 
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(30) 
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where  Pk  =  —u.u.  duijdxj  is  the  turbulence  kinetic  energy  production  term,  cy  =  1.0,  a£  =  1.3,  Cs\ 

=  1.44,  and  Ce 2  =  1.92  are  model  (closure)  constants.  Furthermore,  in  equations  (29)  and  (30),  Gk 
is  the  production  term  relating  to  buoyancy  and  is  given  by 


Gk  =  ~aSi 


vt  dc 
<7  dx. 


(31) 


for  the  case  of  a  two-component  mixture  (where  c  is  the  contaminant  species  concentration).  In 
equation  (30),  Rf  =  - Gk / (Pk  +  Gk)  is  the  flux  Richardson  number  and  C£j  =  0.8  is  an  additional 

model  constant  (Rodi,  1978)  that  is  applicable  both  for  vertical  and  horizontal  buoyant  shear 
layers.  In  particular,  for  a  vertical  shear  layer  where  the  lateral  velocity  component  is  normal  to 
the  direction  of  gravity,  we  note  that  R/=  0.  This  is  the  case  of  relevance  for  our  investigation. 


2.3.2  Validation  test  case:  Thorney  Island 

Extensive  field  trials  on  the  dispersion  of  dense  gas  clouds  at  ground  level  in  the  atmosphere 
were  performed  by  the  Health  and  Safety  Executive  in  Great  Britain.  The  trials  were  conducted 
on  the  Thorney  Island  airfield  in  southern  England  from  July  1982  to  July  1983.  In  these  trials, 
isothermal  gas  clouds  of  fixed  volume  and  of  varying  density  were  instantaneously  released  into 
the  atmosphere.  The  gas  concentrations  in  the  cloud,  consisting  of  Freon- 12  and  air  mixtures, 
were  measured  using  sensors  placed  at  various  downwind  fetches  from  the  release  location.  In 
this  paper,  we  consider  the  trials  conducted  in  the  second  phase  of  the  Thorney  Island  test  series 
(Davies  and  Singh,  1985).  In  these  experiments,  the  dense  gas  cloud  dispersion  was  conducted  in 
the  presence  of  an  obstacle.  The  obstacle  used  was  a  sharp-edged  solid  cube  with  a  dimension  of 
9  m  on  each  side.  In  the  instantaneous  releases,  a  cylindrical  gas  tent  of  14  m  diameter  and  13  m 
height  with  a  total  volume  of  2000  m3  was  filled  with  a  mixture  of  Freon- 12 
(dichlorodifluoromethane,  CCI2F2)  and  nitrogen.  More  specifically,  the  released  gas  consisted  of 
a  mixture  of  68.4%  nitrogen  and  31.6%  Freon-12  (w/w)  to  give  a  mixture  density  that  is  roughly 
twice  that  of  air.  The  top  and  sides  of  the  cylindrical  tent  were  quickly  removed,  leaving  an 
upright  cylinder  of  dense  gas  exposed  to  the  ambient  conditions.  The  mixture  immediately  began 
to  expand  outwards  isotropically  in  all  horizontal  directions  by  gravitational  slumping.  The 
gravitational  slumping  and  atmospheric  dispersion  of  this  dense  gas  cloud  is  affected  by  its 
interaction  with  the  cubical  obstacle,  either  located  downwind  or  upwind  of  the  release. 

In  trial  number  26,  the  windward  face  of  the  cubical  obstacle  was  located  50  m  downwind  from 
the  center  of  the  cylindrical  gas  tent.  The  measured  wind  speed  at  10  m  above  the  ground  was  1.9 
m/s,  and  the  concentration  of  the  gas  cloud  was  measured  on  the  windward  and  leeward  faces  of 
the  obstacle  at  heights  of  6.4  m  and  0.4  m  above  ground  level,  respectively.  In  trial  number  29, 
the  leeward  face  of  the  cubical  obstacle  was  located  27  m  upwind  from  the  center  of  the 
cylindrical  gas  tent.  The  measured  wind  speed  at  10  m  above  the  ground  was  5.6  m/s,  and  the 
concentration  was  measured  on  the  leeward  face  of  the  obstacle  at  the  height  of  0.4  m  above 
ground  level. 

Simulations  of  the  flow  and  dense  gas  dispersion  for  trial  number  26  were  conducted  on  a 
computational  mesh  of  189x60x83  control  volumes  (nodes)  in  a  spatial  domain  (shown  in  Figure 
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1)  with  an  extent  of  16.4//  x  4.4 H  x  14.0//  in  the  streamwise,  vertical  and  spanwise  directions, 
respectively,  where  H  =  9  m  is  the  height  of  the  cubical  obstacle.  The  simulations  for  trial 
number  29  were  performed  on  a  computational  mesh  of  139x60x  83  nodes  in  a  spatial  domain 
that  spanned  15.9 H  x  4 AH  x  14.0 H  in  the  streamwise,  vertical  and  spanwise  directions, 
respectively.  The  mesh  used  for  the  computational  domain  of  trial  number  29  (not  shown)  was 
very  similar  to  that  used  in  trial  number  26,  where  the  grid  lines  were  preferentially  concentrated 
near  the  solid  surfaces  in  order  to  better  capture  the  expected  sharp  gradients  of  the  flow 
properties  here. 

At  the  inlet  boundary  of  the  spatial  domain  used  for  the  simulations,  the  vertical  profile  of  the 
mean  streamwise  velocity  was  approximated  using  a  power-law  functional  form 

u(y)  =  u0  (y  /  y0 )" ,  where  u(j  is  the  reference  velocity  at  the  reference  height  vo  =  10  m.  Table  1 

summarizes  the  values  of  u0  and  n  used  in  the  simulations  for  trial  numbers  26  and  29.  This 

power-law  profile  was  also  used  by  Qin  et  al.  (2007)  and  Sklavounos  and  Rigas  (2004)  for  their 
simulations  of  Thomey  Island.  At  the  outlet  boundary,  zero  Neumann  boundary  conditions  were 
imposed  on  all  flow  variables.  At  the  spanwise  and  top  boundaries,  symmetry  and  free-slip 
boundary  conditions  were  applied  to  the  flow  variables,  respectively.  At  all  the  solid  boundaries 
(ground  and  obstacle  faces),  standard  wall  functions  were  used  for  the  mean  velocities  and 
turbulence  quantities,  and  the  zero  Neumann  boundary  condition  was  applied  for  the 
concentration  of  the  gas  (implying  no  flux  of  material  across  the  solid  boundary). 

In  the  simulations,  computations  were  conducted  firstly  to  predict  the  steady-state  flow  field, 
including  as  such  the  complex  influence  of  both  the  cylindrical  gas  tent  and  the  cubical  obstacle 
on  the  flow.  This  disturbed  flow  field  was  then  used  as  the  initial  condition  for  the  subsequent 
predictions  of  the  dense  gas  dispersion,  which  occurred  immediately  after  the  collapse  of  the 
cylindrical  gas  tent.  The  gas  mixture  inside  the  cylindrical  tent  was  released  instantaneously  to 
the  atmosphere  at  time  t  =  0  s.  Figures  2  and  3  show  the  mean  flow  streamlines  and  contours  of 
concentration  at  different  times  after  the  dense  gas  cloud  was  released  and  documents  the  time 
evolution  of  dense  gas  cloud.  The  gravity  slumping  effect  on  the  dispersion  can  be  seen  clearly 
in  both  figures,  particularly  at  the  earlier  times  where  the  streamlines  within  the  highly 
concentrated  gas  cloud  exhibit  a  strong  downward  bulk  motion.  Note  that  initially  the  buoyancy 
generated  forces  and  pressure  gradients  arising  from  density  differences  between  the  gas  cloud 
and  its  environment  lead  to  a  bulk  motion  that  causes  the  dense  cloud  to  spread  in  all  directions 
near  the  release  location  (including  a  lateral  spreading,  as  well  as  an  upwind  spreading  against 
the  prevailing  wind  direction).  Further  downstream  from  the  release,  the  bulk  flow  resulting  from 
the  gravitational  slumping  weakens  as  diffusion,  mixing,  and  entrainment  between  the  gas  cloud 
and  the  ambient  air  reduces  the  negative  buoyancy  effects  of  the  cloud  and  strengthens  the 
influence  of  the  externally  imposed  velocity  field  on  the  transport  and  dispersion  of  the  cloud.  It 
is  noted  in  trial  number  26  that  the  dense  gas  cloud  flows  around  and  over  the  cubical  obstacle 
located  downwind  of  the  release  (cf.  Figure  2).  In  contrast,  for  trial  number  29,  the  cubical 
obstacle  located  upwind  of  the  release  is  seen  to  block  the  further  upwind  spread  of  the  dense  gas 
cloud  arising  from  the  gravitational  slumping.  Also,  in  addition  to  the  blocking  effect,  the  gas 
cloud  undergoes  increased  dilution  in  the  wake  of  the  upstream  cubical  obstacle  (cf.  Figure  3),  as 
well  as  retention  or  “hold-up”  of  the  material  within  the  wake  of  the  obstacle  followed  by  a  slow 
detrainment  of  material  from  the  wake.  This  retention  of  material  in  the  wake  followed 
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subsequently  by  the  detrainment  of  this  material  converts  the  putative  instantaneous  release  into 
a  time- varying  one. 

Figures  4  and  5  display  a  comparison  between  the  measured  and  predicted  gas  concentration 
time  histories  for  trial  numbers  26  and  29,  respectively.  In  particular,  Figure  4  shows  that  the 
time  histories  of  concentration  observed  at  the  windward  and  leeward  faces  of  the  cubical 
obstacle  are  quite  different.  Note  that  the  persistence  of  the  gas  cloud  (as  measured  by  the 
difference  between  the  arrival  and  departure  times  of  the  cloud)  observed  at  the  leeward  face  of 
the  cubical  obstacle  is  significantly  greater  than  that  observed  at  the  windward  face.  This 
demonstrates  the  effect  of  the  obstacle  on  dispersion,  where  the  recirculation  bubble  in  the  wake 
region  of  the  obstacle  entrains  and  “holds-up”  (traps)  the  contaminant  material  of  the  gas  cloud. 
The  hold-up  mechanism  within  the  wake  of  the  obstacle  is  correctly  predicted  by  the  model. 

Generally,  numerical  predictions  of  concentration  are  seen  to  be  in  quite  good  agreement  with 
the  experimental  data.  It  is  noted  that  in  interpreting  the  level  of  agreement  between  predictions 
and  measurements,  one  needs  to  recognize  that  the  measurements  of  the  short  time-averaged 
concentration  obtained  in  trial  numbers  26  and  29  correspond  to  one  realization  of  the 
instantaneous  dispersion,  whereas  the  predictions  correspond  to  an  ensemble-averaged 
concentration  (and,  as  such,  is  associated  with  an  average  over  an  ensemble  of  realizations  of  the 
instantaneous  dispersion).  Because  the  model  prediction  is  compared  to  one  realization  out  of  the 
ensemble  of  all  possible  realizations,  the  model  cannot  be  expected  to  provide  greater  precision 
in  its  predictions  than  the  expected  variability  in  the  atmospheric  dispersion.  Keeping  this  caveat 
in  mind,  the  discrepancies  between  the  predictions  and  the  measurements  could  be  attributed  to 
the  use  of  the  standard  k-s  turbulence  model.  It  is  well  known  that  the  standard  k-s  model 
tends  to  over-predict  levels  of  the  turbulence  kinetic  energy  for  a  flow  impinging  on  a  solid 
obstacle  [the  so-called  the  “stagnation  point  anomaly”,  (Durbin,  1996)],  resulting  in  prediction  of 
a  larger  recirculation  bubble  in  the  wake  region  of  a  three-dimensional  obstacle  than  is  actually 
observed. 

For  trial  number  26,  the  more  pronounced  impingement  of  the  approaching  flow  on  the  front  face 
of  the  cylindrical  gas  tent  resulting  in  a  predicted  larger  recirculation  bubble  in  the  wake  region 
of  the  tent  will  have  the  tendency  to  dilute  the  gas  concentration  faster  after  the  tent  collapses.  In 
consequence,  this  will  lead  to  a  lower  predicted  concentration  downstream  of  the  cylindrical  gas 
tent  when  compared  to  the  associated  measured  concentration,  as  can  be  seen  in  Figure  4.  For 
trial  number  29,  a  larger  recirculation  bubble  behind  the  upwind  cubical  obstacle  will  entrain 
material  from  the  dense  gas  cloud  (that  has  spread  upwind  owing  to  the  gravity  slumping)  into 
the  wake  earlier.  This  effect  will  result  in  an  earlier  predicted  arrival  time  for  the  cloud  (at  the 
leeward  face  of  the  upwind  obstacle)  in  comparison  to  the  measurements,  as  can  be  seen  in 
Figure  5. 

The  effect  of  the  buoyancy  production  term  Gk  in  the  k-s  model  [cf.  equations  (29)  and  (30)] 
on  the  prediction  of  dense  gas  dispersion  is  investigated  and  demonstrated  in  Figure  6.  It  was 
found  that  the  inclusion  of  Gk  term  in  the  k-s  model  suppresses  the  turbulence  kinetic  energy, 
particularly  in  the  wake  region  of  the  obstacle,  which  is  not  too  surprising  owing  to  the  negative 
buoyancy  effect  of  the  dense  gas  cloud.  In  contrast,  the  effect  of  the  Gk  term  on  the  prediction  of 
the  concentration  levels  in  the  cloud  appears  to  be  insignificant.  This  observation  is  similar  to  the 
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findings  of  Worthy  et  al.  (2001),  where  they  reported  that  the  Gt  term  made  very  little  difference 
to  the  modeling  of  buoyant  thermal  plumes. 


3.  DESCRIPTION  OF  SCENARIOS  -  DOWNWIND  HAZARD  MODELING 

The  downwind  hazard  modeling  for  the  multi-spectrum  scenario  will  consist  of  the  following 
components,  and  will  involve  two  primary  locations  in  Calgary;  namely  along  the  railway  tracks 
south  of  the  intersection  of  (1)  9th  Avenue  SE  and  4th  Street  SE  and  (2)  9th  Avenue  SW  and  11th 
Street  SW.  The  release  involving  the  van  will  take  place  within  about  1  block  of  the  “accidental” 
release  along  the  railway  line.  Suggested  locations  relative  to  the  primary  release  along  the 
railway  line  are  summarized  below.  Note  that  the  release  involving  chlorine  and  a  chemical 
agent  each  involves  two  components  -  an  near-instantaneous  release  and  a  longer  continuous 
release  involving  evaporation  of  the  material.  During  the  evaporation  phase,  first  responders  with 
the  proper  detection  equipment  may  be  able  to  identify  the  chemical  agent  released  from  the 
nearby  van.  The  toxic  industrial  materials  can  be  identified  probably  from  UN  labels  on  the 
railway  tanker  cars. 


3.1  Scenario  1:  derailment  of  a  rail  car  containing  methanol  (UN  1230,  Class  3)  -  liquid 
density  of  methanol  is  791  kg/m3 

Rail  car  could  potentially  release  up  to  34,500  US  gallons  of  methanol.  It  will  be  assumed  that 
the  rail  car  is  only  %  full,  so  the  rail  car  only  holds  about  25,000  US  gallons  of  methanol. 
Methanol  is  flammable,  and  it  is  assumed  that  the  process  causing  the  derailment  also  results  in  a 
fire.  The  fire  is  assumed  to  destroy  about  20%  of  the  methanol.  Of  the  remaining  methanol  not 
destroyed  by  the  fire  (about  20,000  US  gallons,  or  about  75  m3),  roughly  A  of  this  volume  was 
“flashed”  to  the  atmosphere  in  form  of  vapor/aerosol  droplets  and  remaining  methanol  was 
released  to  atmosphere  in  continuous  manner  at  a  constant  emission  rate.  Initial  “flash”  caused 
instantaneous  release  of  791x37.5  ~  30,000  kg  of  methanol.  The  remaining  methanol  was 
“boiled”  off  by  the  fire  at  a  constant  emission  rate  of  1000  kg/min  (for  about  30  min). 

Location  1:  51  deg,  02’,  37.12”  N,  114  deg,  03’,  11.89”  W 
Location  2:  51  deg,  02’,  41.07”  N,  114  deg,  05’,  14.09”  W 

Note  -  Location  1  is  along  the  rail  tracks  that  lie  south  of  9th  Avenue  SE.  The  location  lies  about 
1  block  south  of  the  intersection  of  9th  Avenue  SE  and  4th  Street  SE.  Location  2  is  along  the  rail 
tracks  that  lie  south  of  9th  Avenue  SW.  The  location  lies  about  1  block  south  of  the  intersection 
of  9th  Avenue  SW  and  1 1th  Street  SW,  along  the  rail  tracks  that  run  through  downtown  Calgary. 


3.2  Scenario  2:  derailment  of  a  rail  car  containing  chlorine  (UN  1017,  Class  2) 

Rail  car  (model  ERG  ID  1017)  which  can  hold  up  to  25,700  US  gallons  of  liquefied  chlorine  is 
only  %  full.  It  is  assumed  that  the  derailment  caused  10,000  kg  of  chlorine  to  be  expelled  into  the 
atmosphere  in  the  form  of  vapor/aerosol  droplets.  This  material  was  “flashed”  to  the  atmosphere 
over  a  short  time  of  100  s,  at  a  release  rate  of  100  kg/s  over  this  period  of  time.  In  addition,  5,000 
kg  of  liquid  chlorine  was  spilled  on  the  ground  as  a  result  of  the  derailment.  The  liquid  pool  of 
chlorine  is  assumed  to  be  evaporating  from  the  ground  surface  at  a  constant  rate  of  2.78  kg/s  over 
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a  period  of  30  minutes,  at  which  time  all  chlorine  in  the  liquid  pool  would  have  evaporated.  Of 
course,  there  is  still  quite  a  bit  of  liquid  chlorine  that  remains  inside  the  rail  car  that  will  need  to 
be  cleaned  up  by  the  first  responders. 

Location  1:  51  deg,  02’,  37.12”  N,  114  deg,  03’,  11.89”  W 
Location  2:  51  deg,  02’,  41.07”  N,  114  deg,  05’,  14.09”  W 

Locations  1  and  2  are  the  same  as  those  in  Scenario  1 . 

For  the  simulations  of  Scenarios  1  and  2  above,  several  wind  directions,  including  45°,  90°, 
135°,  225°,  270°  and  315°  are  chosen  that  will  carry  the  material  into  the  downtown  area  of 
Calgary  (where  all  the  tall  office  towers  are). 

3.3  Scenario  3:  van  containing  two  drums  of  an  unknown  chemical  agent 

'J 

This  scenario  corresponds  to  a  volume  of  about  0.38  m  .  Terrorists  have  planted  explosives  in 
the  van,  and  when  the  explosive  incident  occurs,  it  is  assumed  about  a  certain  quantity  of  the 
chemical  agent  is  destroyed  in  the  explosive  event.  A  specified  amount  of  the  chemical  agent  is 
released  instantaneously  as  a  result  of  the  explosion  in  the  form  of  vapor/aerosol  droplets,  and 
the  remaining  amount  of  liquid  chemical  agent  spilled  on  the  ground.  The  material  evaporates 
from  the  ground  surface  at  a  constant  rate  of  over  a  period  of  time.  The  actual  quantities  released 
are  not  reported  here.  They  are  available  in  a  separate  confidential  document  (private 
communication  Yee,  2012). 

The  location  of  the  van  is  with  a  block  or  two  of  the  site  of  the  derailment,  and  9  possible 
locations  are  chosen  as  follows: 

Relative  to  “accident”  at  Location  1,  the  following  4  locations  are  used  in  our  simulations,  in 
which  only  the  southeast  (or  135°)  wind  direction  is  considered: 

1.  51  deg,  02’,  39.18”  N,  114  deg,  03’,  17.20”  W  (referred  to  as  L3-1) 

2.  51  deg,  02’,  38.62”  N,  114  deg,  03’,  20.18”  W  (referred  to  as  L3-2) 

3.  51  deg,  02’,  35.87”  N,  114  deg,  03’,  16.91”  W  (referred  to  as  L3-3) 

4.  51  deg,  02’,  35.81”  N,  114  deg,  03’,  19.56”  W  (referred  to  as  L3-4) 

Relative  to  “accident”  at  Location  2,  the  following  5  locations  are  used  in  our  simulations,  in 
which  only  the  southwest  (or  225°)  wind  direction  is  considered: 

1.  51  deg,  02’,  41.46”  N,  114  deg,  05’,  21.06”  W  (referred  to  as  L4-1) 

2.  51  deg,  02’,  41.30”  N,  114  deg,  05’,  17.71”  W  (referred  to  as  L4-2) 

3.  51  deg,  02’,  40.57”  N,  1 14  deg,  05’,  20.43”  W  (referred  to  as  L4-3) 

4.  51  deg,  02’,  40.50”  N,  1 14  deg,  05’,  18.24”  W  (referred  to  as  L4-4) 

5.  51  deg,  02’,  40.96”  N,  114  deg,  05’,  19.48”  W  [referred  as  L4-LC  (level  crossing 
plume  point)]. 
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4.  RESULTS  AND  DISCUSSION  FOR  DIFFERENT  SCENARIOS 


The  modeling  (outer)  domain  with  the  extent  of  6,400  m  x  4,200  m  x  830  m  in  the  x-  (or,  W-E), 
y-  (or,  S-N)  and  z-  (or,  vertical)  directions,  respectively,  with  its  southwest  (SW)  and  northeast 
(NE)  comers  located  at  (702,260  m  E,  5,656,815  m  N)  and  (708,660  m  E,  5,661,015  m  N), 
respectively,  in  the  Universal  Transverse  Mercator  (UTM)  coordinate  system  are  shown  in 
Figures  7  and  8.  The  internal  (dimensionless)  coordinate  system  used  in  urbanSTREAM  is 
normalized  by  a  reference  length  scale  which  is  chosen  in  this  case  to  be  L=  3,593  m.  A  proper 
sub-domain  within  this  modeling  region  is  chosen  as  the  “inner  region”  in  which  buildings  will 
be  explicitly  resolved  in  the  flow  simulation.  For  this  example,  this  rectangular  building-aware 
region  of  the  size  of  3,600  m  x  2,100  m  has  its  southwest  comer  at  (703,660  m  E,  5,657,865  m 
N)  and  its  northeast  comer  at  (707,260  m  E,  5,659,965  m  N)  in  the  UTM  coordinate  system.  In 
the  portion  of  the  solution  domain  lying  outside  the  building-aware  (inner)  region,  all  buildings 
are  treated  as  virtual  and  their  effects  on  the  flow  are  modeled  using  a  distributed  drag  force 
representation  in  the  mean  momentum  equations  [see,  e.g.,  Lien  &  Yee  (2004,  2005),  Lien  et  al. 
(2005)]. 

A  building  database  in  the  form  of  an  Arc  View  Shapefile  for  the  City  of  Calgary  is  initially 
available  from  Environment  Canada.  However,  it  was  identified  in  the  early  stage  of  the  project 
that  a  few  key  buildings,  such  as  the  Bow  Building  (see  Figure  9),  was  missing  in  the  database. 
A  new  building  database  was  later  on  provided  by  3DIntemet  Inc.  and  used  for  the  rest  of  our 
CFD  simulations.  These  Shapefiles  were  used  in  urbanGRID  to  generate  automatically  a  grid 
mesh  over  the  modeling  region  as  shown  in  Figure  10.  Then,  the  simulation  was  carried  out  in  a 
3-D  Cartesian  framework,  and  curved  surfaces  on  buildings  or  planar  building  surfaces  that  are 
not  aligned  with  the  grid  lines  are  approximated  by  stepwise  surfaces.  A  mesh  of  498x336x65 
grid  lines  in  the  x-,  y-,  and  z-directions,  respectively,  was  used  to  accommodate  all  the  necessary 
geometrical  details.  The  inner  building-aware  region  was  covered  with  a  fine  mesh  of 
450x300x65  grid  lines  to  better  approximate  the  building  features  in  this  region.  The  grid 
arrangement  adopted  here  is  shown  on  the  x-y  plane  in  Figure  8.  Hence,  the  fine  grid  used  in  the 
building-aware  region  contains  8,775,000  nodes,  whereas  the  entire  computational  domain  was 
covered  with  a  mesh  of  10,876,320  nodes.  The  grid  lines  were  preferentially  concentrated  near 
the  solid  surfaces  (ground,  building  rooftops  and  walls)  where  the  gradients  in  the  flow 
properties  are  expected  to  be  the  greatest,  and  the  spacings  between  the  grid  lines  were  gently 
stretched  with  an  increasing  distance  from  the  solid  surfaces. 

The  flow  field  in  the  computational  domain  was  computed  using  urbanSTREAM  in  a  standalone 
mode,  and  the  inflow  condition  for  a  given  wind  direction  was  described  in  Section  2.2.  The  flow 
simulation  in  the  computational  domain  was  carried  out  in  the  steady  RANS  mode  assuming 
there  was  no  change  in  the  mean  speed  or  mean  direction  of  the  flow  within  the  45  minutes 
simulation  time  period.  The  effects  of  the  buildings  on  the  flow  field  outside  the  building-aware 
region  (viz.,  in  the  “virtual  building”  region)  were  represented  by  a  distributed  drag  force 
approximation  in  the  mean  momentum  equation  given  below: 
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with  a  normalized  drag  coefficient  of  Cd  =  100.  The  specification  of  this  value  is  based  on  our 
earlier  calculations  for  the  IOP9  test  problem  in  comparison  with  the  experimental  data  [see,  e.g., 

Lien  et  al  (2008)].  In  equation  (33),  A  is  the  frontal  area  density  (frontal  area  of  buildings 
exposed  to  the  wind  per  unit  volume)  and  Q  =  \u\  is  the  magnitude  of  the  spatially-averaged 

time-mean  wind  speed.  A  modified  k  —  s  model  consistent  with  equations  (32)-(33)  can  be 
found  in  Lien  et  al.  (2005).  Sample  results  at  the  wind  directions  of  135°  and  225°  within  and  in 
the  building-resolved  regions  of  the  city  are  shown  in  Figures  11  and  12  to  illustrate  the 
complexity  of  streamline  patterns  and  associated  velocity  fields.  Most  results  presented  in 
Section  4  are  obtained  with  the  parallel  version  of  urbanSTREAM  (or  urbanSTREAM-P) 
running  on  SHARCNET  using  16  nodes  unless  stated  otherwise. 


4.1  Scenario  1:  derailment  of  a  rail  car  containing  methanol 

In  Scenario  1,  the  amount  of  methanol  in  the  rail  car  is  34,500  galx(3/4)  [3/4  full]=25,875  gal,  in 
which  only  80%  of  methanol  is  left  [20%  of  the  methanol  is  destroyed];  i.e.,  only  25,875  gal 
x4/5  =20,700  gal  (=78.35  m3)  of  methanol  is  left.  Of  the  remaining  methanol  not  destroyed  by 
the  fire,  we  can  split  the  dispersion  processes  into  two  parts: 

1.  “Instantaneous  puff  release”  (“puff  release”  or  R2  release):  Vi  of  this  volume  [i.e., 
20,700x1/2=10,350  gal  (or  39.18  m3)]  was  “flashed”  to  the  atmosphere  in  form  of 
vapor/aerosol  droplets.  The  density  of  methanol  at  standard  temperature  and  temperature 
(STP)  condition  is  about  791  kg/  m3.  So  39.18x791=30,991  kg  ~  30,000  kg  of  methanol, 
which  is  treated  as  a  neutrally  buoyant  (or  passive)  gas  in  our  simulations,  was  “flashed”  to 
the  atmosphere  during  the  first  time  step:  At  =  10  sec. 

2.  “Continuous  release”  (or  R1  release):  another  Vi  of  this  volume  “boiled”  off  by  the  fire  for 
about  30  min  at  a  constant  emission  rate  of  30,000/30=1000  kg/min=16.667  kg/s. 

Since  both  releases  above  are  governed  by  the  linear  advection-diffusion  equation  for 
concentration  of  methanol  with  a  given  velocity  field,  the  superposition  principle  applies. 
Therefore,  simulations  for  continuous  and  puff  releases  are  done  separately,  and  the  resulting 
concentration  fields  to  be  shown  later  are  the  solutions  obtained  by  combining  both  releases 
together. 

Sample  streamline  patterns  superimposed  with  iso-surfaces  of  concentration  at  225°  wind 
direction  are  shown  in  Figures  13  (front  view)  and  14  (rear  view),  respectively,  to  illustrate  the 
complexity  of  flow  and  concentration  fields  in  the  wake  region  of  buildings. 
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Figures  15  to  20  are  contours  of  concentration  on  a  logarithmic  scale  for  the  source  at  Location  1 
(51  deg,  02’,  37.12”  N,  114  deg,  03’,  11.89”  W)  at  t  =  5  and  30  min  for  the  wind  direction  of 
45°,  90°  and  135°,  respectively.  For  the  source  at  Location  2  (51  deg,  02’,  41.07”  N,  114  deg, 
05’,  14.09”  W),  similar  contours  of  concentration  are  given  in  Figures  21  to  26  at  t  =  5  and  30 
min  for  the  wind  direction  of  225°,  270°  and  315°.  In  these  figures,  iso-surfaces  of  logio(C)  =  - 
2.2,  -3.3,  -4.4,  where  the  unit  of  C  is  kg/m3,  are  shown,  and  logio(C)  <  12  are  omitted  because 
they  are  insignificant  for  practical  applications. 

Let  us  consider,  e.g.,  Figures  15  and  16  at  t  =  5  and  30  min  for  an  incident  wind  direction  of  45°. 
It  is  clearly  seen  that  both  R1  and  R2  releases  co-exist  in  Figure  15  in  terms  of  logio(C)=-4.4  iso¬ 
surface.  However,  at  t  =  30  the  “puff’  (associated  with  the  R2  release)  already  leaves  the  solution 
domain,  and  only  iso-surfaces  for  R1  (continuous)  release  are  present  in  Figure  16.  Similar 
comments  apply  to  the  figures  at  the  other  wind  directions  (i.e.,  Figures  17  to  20  with  a  source  at 
Location  1,  and  Figures  21  to  26  with  a  source  at  Location  2). 

Before  discussing  results  for  Scenario  2,  we  repeat  the  calculations  for  Scenario  1  except 
replacing  methanol  by  liquefied  chlorine,  and  using  the  heavy-gas  model  described  in  Section 
2.3  to  simulate  dispersion  of  chlorine  with  the  source  at  Location  1  at  the  wind  direction  of  135°. 

Let  us  first  introduce  the  following  dimensionless  groups: 
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where  the  reference  length  L  =  3,593  m,  and  the  reference  velocity  U  =5.556  m/s  in  the  present 
calculations.  Equation  (14)  in  dimensionless  form  can  be  written  as 
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where  Re  =  P°^  ,  and 
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In  equation  (36),  M0  =28.97  kg/kmol  and  Mg  =70.91  kg/kmol  are  the  molecular  weights  for 

air  and  chlorine,  respectively.  Similarly,  equation  (15)  combined  with  equation  (28)  can  be  non- 
dimensionalized  as 
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where  Q  in  kg/s  is  the  source  strength  and  xs  is  the  source  location.  Note  that  the  source  term  in 

equation  (37)  was  not  in  equation  (15),  and  the  unit  for  S(x  —  xs)  is  m  3.  To  be  consistent  with 

the  Thomey  Island  test  problem,  crc  =  0.63  is  chosen  here,  which  is  the  turbulent  Schmidt 
number.  Also,  the  molecular  diffusion  in  equation  (37)  is  ignored  in  order  to  be  compatible  with 
the  standard  (or  high-Re)  k  —  s  model  adopted  in  the  present  study.  Integrating  the  last  term  of 
equation  (37)  over  a  (dimensionless)  finite  volume  containing  the  source  on  the  ground  gives  us 
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where  Q  =  16.667  kg/s .  cs  in  equations  (34)  and  (38)  is  chosen  to  be  1  kg/m  for  simplicity.  In 
the  following,  two  approaches  will  be  introduced  into  urbanEU  to  model  the  (effect  of)  source 
term:  one  through  Q  and  the  other  through  the  concentration  boundary  condition  on  the  ground 
at  x  =  xs . 


4.1.1  Q-approach 


In  the  “Q-approach”, 


Q 

cfiU 


in  urbanEU  is  turned  on  when  0  <  t  <  30  min.  The  simulation  starts 


with  the  use  of  urbanSTREAM  to  obtain  a  quasi-steady  state  velocity  and  turbulence  field  for  a 
constant-density  (incompressible)  flow  (viz.,  the  solution  converges  at  around  10  iterations  per 
time  step).  With  this  as  the  initial  condition,  urbanEU  is  then  coupled  to  urbanSTREAM  to 
restart  the  calculations  for  another  270  time  steps  at  At  =  10  sec. 


4.1.2  Wjet-approach 


To  further  validate  the  implementation  of  heavy- gas  modeling  capability  in  urbanEU, 
alternative  approach  -  referred  to  as  the  “Wjet  approach”  is  proposed  here.  Firstly,  we 

Q 


cL2U 


=  0 .  At  the  source  location  (i.e.,  at  x  =  xs  )  on  the  ground,  we  set 


an 

set 


1 .  c  =  1  (or  c  =  cs ), 


2. 


w  = 


Q 

csUAs  ’ 


where  w  is  the  dimensionless  vertical  component  of  the  (jet)  velocity  at  the  source  location,  and 
As  is  the  source  surface  area. 


Sample  results  showing  iso-surfaces  of  concentration  of  log(C)  =  -5,  -7  and  -10  at  t  =  10  min 
obtained  with  both  the  Q-  and  Wjet  approaches  are  shown  in  Figures  27  and  28.  It  is  observed 
that  both  results  are  virtually  identical. 
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We  also  compare  contours  of  log(C)  superimposed  with  the  velocity  vector  field  near  the  source 
location  at  t  =  10  min  obtained  with  the  dense-gas  model  (using  the  Q-approach)  and  the  passive- 
gas  model,  which  are  shown  in  Figures  29  and  30,  respectively.  Interestingly,  both  results  are 
quite  similar.  It  is  noted  here  that  in  the  case  of  the  passive-gas  approach,  urbanEU  is  decoupled 
from  urbanSTREAM,  and  only  the  concentration  equation  is  required  to  be  solved.  In  contrast,  if 
the  dense-gas  model  is  employed,  urbanEU  and  urbanSTREAM  are  strongly  coupled  together, 
and  typically  175  “outer”  iterations  per  time  step  is  required  to  reach  convergence.  In  terms  of 
CPU  time  consumption,  the  passive-gas  approach  is  approximately  400  times  faster  than  the 
dense-gas  approach  (whether  using  either  the  Q-approach  or  Wjet  approach).  Since  the  results 
shown  in  Figures  29  and  30  seem  to  suggest  that  there  is  no  major  difference  using  either  the 
dense-gas  or  the  passive-gas  approach,  only  the  passive-gas  approach  will  be  considered  for  the 
rest  of  the  scenarios  to  save  computing  resources  as  there  are  many  cases  that  need  to  be 
computed  for  the  present  project. 

4.2  Scenario  2:  derailment  of  a  rail  car  containing  chlorine 

In  Scenario  2,  the  amount  liquefied  chlorine  of  in  the  rail  car  (model  ERG  ID  1017)  is  25,700 
galx(l/2)  [1/2  full]=12,850  gal  ~  76,000  kg  based  on  the  density  of  the  liquefied  chlorine  being 
1562.5  kg/m3  (at  1.013  bar  at  boiling  point).  It  is  assumed  that  the  derailment  caused  10,000  kg 
of  chlorine  to  be  expelled  into  the  atmosphere  in  the  form  of  vapor/aerosol  droplets  over  a  short 
time  of  100  sec,  at  a  release  rate  of  10,000/10=100  kg/s  over  this  period  of  time,  which  will  be 
considered  as  “instantaneous  release”  (or  R2  release)  in  our  simulations. 

In  addition,  5,000  kg  of  liquid  chlorine  that  was  spilled  on  the  ground  is  assumed  to  be 
evaporating  from  the  ground  surface  over  a  period  of  30  minutes  (or  1,800  sec)  at  a  constant  rate 
of  5,000/1,800  =2.78  kg/s,  at  which  time  all  chlorine  in  the  liquid  pool  would  have  evaporated, 
and  treated  as  “continuous  release”  (or  R1  release)  in  our  simulations.  In  both  releases,  chlorine 
gas  is  considered  as  a  passive  gas  to  save  computing  resources. 

Figures  31  to  36  are  contours  of  concentration  on  a  logarithmic  scale  for  the  source  at  Location  1 
(51  deg,  02’,  37.12”  N,  114  deg,  03’,  11.89”  W)  at  t  =  5  and  30  min  for  the  wind  direction  of 
45°,  90°  and  135°.  For  the  source  at  Location  2  (51  deg,  02’,  41.07”  N,  114  deg,  05’,  14.09”  W), 
similar  contours  of  concentration  are  given  in  Figures  37  to  42  at  t  =  5  and  30  min  for  the  wind 
direction  of  225°,  270°  and  315°.  In  these  figures,  iso-surfaces  of  logio(C)  =  -2.2,  -3.3,  -4.4  are 
shown,  and  logio(C)  <  12  are  omitted  for  the  same  reason  as  in  Scenario  1. 

Note  that  the  3  iso-surfaces  in  Figures  15  to  26  for  Scenario  1,  and  in  Figures  31  to  42  for 
Scenario  2  are  chosen  at  the  same  values  of  logio(C)  =  -2.2,  -3.3,  -4.4  for  comparison  purposes 
later.  Let  us  compare  R1  (continuous)  and  R2  (instantaneous  or  puff)  releases  in  Scenarios  1  and 
2  based  on  the  following  Table: 


R1  release 

R2  release 

Scenario  1 

30,000  kg  over  30  min 

30,000  kg  over  10  sec 

Scenario  2 

5,000  kg  over  30  min 

10,000  kg  over  100  sec 
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Similar  to  what  we  observe  in  the  results  for  Scenario  1,  iso-surfaces  of  logio(C)  for  both  R1  and 
R2  releases  of  Scenario  2  co-exist  at  t  =  5  min  for  all  wind  directions  with  a  source  at  Location  1 
as  shown  in  Figures  31,  33  and  35,  and  with  a  source  at  Location  2  as  shown  in  Figures  37,  39 
and  41.  In  contrast,  at  t  =  30  min  the  “puff’  (or  R2  release)  has  been  convected  out  of  the 
solution  domain,  and  only  iso-surfaces  associated  with  the  R1  release  can  be  seen  in  Figures  32, 
34,  36,  38,  40  and  42.  Also,  at  t  =  30  min  for  a  given  wind  direction,  say  135°,  the  volume 
corresponding  to  logio(C)  =  -4.4  for  Scenario  2  (shown  in  Figure  36)  is  smaller  than  that  for 
Scenario  1  (shown  in  Figure  20).  This  is  because  the  amount  of  R1  release  (30,000  kg)  for 
Scenario  1  is  6  times  larger  than  that  (5,000  kg)  for  Scenario  2. 


4.3  Scenario  3:  van  containing  two  drums  of  an  unknown  chemical  agent 

In  Scenario  3,  it  is  assumed  that  a  quantity  of  an  unknown  chemical  agent  is  released 
instantaneously  as  a  result  of  the  explosion  in  the  form  of  vapor/aerosol  droplets.  This  is  treated 
as  “instantaneous  puff  release”  during  the  first  time  step  of  At  =  10  sec. 

In  addition,  as  a  result  of  the  explosion  in  the  form  of  vapor/aerosol  droplets,  and  a  certain 
quantity  of  the  liquid  agent  is  spilled  on  the  ground.  The  material  evaporates  from  the  ground 
surface  at  a  constant  rate  over  a  period  of  30  minutes. 

Four  source  locations  close  to  “Location  1”  (at  51  deg,  02’,  37.12”  N,  114  deg,  03’,  11.89”  W) 
for  both  Scenarios  1  and  2,  referred  to  as  L3-1,  L3-2,  L3-3  and  L3-4  in  Section  3.3,  are  shown  in 
Figures  43,  47,  51  and  55  on  a  google  map.  The  corresponding  top  (or  x-y)  view  of  contours  of 
logio(C),  where  C  is  concentration  in  arbitrary  units,  are  given  in  Figures  44,  48,  52  and  56  at  a 
wind  direction  of  135°  (or  southeast  wind).  Figures  45  and  46,  Figures  49  and  50,  Figures  53  and 
54,  and  Figures  57  and  58  are  contours  of  logio(C)  in  conjunction  with  3  iso-surfaces  of 
logio(C)=-5,  -6  and  -7  (with  C  in  arbitrary  units)  at  t  =  5  and  15  min  for  sources  at  L3-1,  L3-2, 
L3-3  and  L3-4,  respectively. 

Similarly,  five  source  locations  close  to  “Location  2”  (at  51  deg,  02’,  41.07”  N,  114  deg,  05’, 
14.09”  W)  for  Scenarios  1  and  2,  referred  to  as  L4-1,  L4-2,  L4-3,  L4-4  and  L4-LC  in  Section  3.3, 
are  shown  in  Figures  59,  63,  67,  71  and  75  on  a  google  map.  The  corresponding  top  (or  x-y)  view 
of  contours  of  logio(C)  are  given  in  Figures  60,  64,  68,  72  and  76  at  a  wind  direction  of  225°  (or 
southwest  wind).  Figures  61  and  62,  Figures  65  and  66,  Figures  69  and  70,  Figures  73  and  74, 
and  Figures  77  and  78  show  contours  of  logio(C)  in  conjunction  with  3  iso-surfaces  of  logio(C)=- 
5,  -6  and  -7  at  t  =  5  and  15  min  for  sources  at  L4-1,  L4-2,  L4-3,  L4-4  and  L4-LC,  respectively. 

As  seen  from  Figures  44,  48,  52  and  56  for  sources  located  close  to  Location  1  (at  51  deg,  02’, 
37.12”  N,  114  deg,  03’,  11.89”  W)  for  a  southeast  wind  (i.e.,  at  a  wind  direction  of  135°),  the 
“domain  of  influence”  [or  spread  of  the  iso-surfaces  of  logio(C)]  downstream  of  the  source  is 
mainly  determined  by  the  wind  direction  and  source  location.  However,  the  concentration  fields 
near  sources  at  L3-1  and  L3-4  (see  Figures  44  and  56)  are  strongly  influenced  by  the  presence  of 
buildings  immediately  upstream  of  the  sources.  Let  us  examine  contours  and  iso-surfaces  of 
logio(C)  at  t  =  5  and  15  min  shown  in  Figures  45  and  46  at  L3-1  source  location.  It  is  seen  that 
both  R1  and  R2  releases  still  co-exist  in  the  solution  domain  even  at  t  =  15  min.  However, 
contours  projected  on  the  y-z  and  x-z  planes  suggest  that  the  “puff’  is  moving  (and  decaying)  in 
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the  southeast  direction,  which  is  consistent  with  the  wind  direction  of  135°.  Similar  observation 
applies  to  L3-2,  L3-3  and  L3-4  source  locations. 

For  sources  located  close  to  Location  2  (51  deg,  02’,  41.07”  N,  114  deg,  05’,  14.09”  W)  at  a 
wind  direction  of  225°,  the  top  view  of  the  overall  shape  of  the  downstream  plume  in  terms  of 
logio(C)=-7  is  quite  similar  for  all  five  source  locations  examined  here  (see  Figures  60,  64,  68, 
72  and  76).  This  suggests  that  as  the  size  of  the  plume  grows  to  the  extent  that  its  length  scale  is 
larger  than  the  average  building  height,  the  shape  of  the  plume  is  hardly  affected  by  the  detailed 
configuration  (i.e.,  shape,  size  and  orientation)  of  the  buildings  underneath.  The  basic  idea  of  the 
drag-force  model  described  by  equations  (32)  and  (33)  is  essentially  based  on  the  above 
observation,  which  hypothesizes  that  the  effects  of  buildings  only  influence  the  large  scale  of  the 
wind  field  and  concentration  fields  [after  volume-averaging  of  the  RANS  equation  (see,  e.g., 
Lien  et  al.,  2005)].  Note  that  all  buildings  in  the  “inner  region”  are  explicitly  resolved  here.  Also, 
contours  projected  on  the  y-z  and  x-z  planes  in  Figures  62,  66,  70,  74  and  78  at  L4-1,  L4-2,  L4-3, 
L4-4  and  L4-LC  source  locations  suggest  that  the  “puff’  is  moving  (and  decaying)  in  the 
southwest  direction,  which  is  consistent  with  the  wind  direction  of  225°. 


4.4  Development  of  a  compact  data  structure  for  file  transfer 

To  facilitate  the  transfer  of  concentration  data  over  40  min  at  the  time  interval  of  every  1  sec  for 
0<  t<l  min,  every  5  sec  for  1  min  <  t  <10  min,  and  every  10  sec  for  10  min  <  t  <40  min  on  a 
grid  of  the  size  of  ~  11  million  nodes,  a  compact  data  structure  is  developed  for  this  project,  in 
which  only  the  concentration  values  of  greater  than  the  cut-off  number  (chosen  to  be  10~20  kg/m3 
in  the  present  report)  are  stored. 


4.4.1  Data  structure  and  its  parallel  implementation 

In  urbanSTREAM-P,  a  new  module:  datatype .  f  90  is  developed  and  listed  below: 

MODULE  datatype 
!  »TYPE 


TYPE  OUTPUT  TYPE 

INTEGER 

:  MX 

INTEGER 

:  IL,  JL , 

KL 

INTEGER 

:  IG,  JG, 

KG 

END  TYPE 

OUTPUT  TYPE 

!  »TYPE 
!  »TYPE 

TYPE  OUTPUTC_TYPE 

INTEGER  ::  I,  J,  K 
REAL : :  C 

END  TYPE  OUTPUTC_TYPE 
!  »TYPE 

TYPE (OUTPUT_TYPE  ),  allocatable,  dimension )  ::  OUTPUT 

TYPE (OUTPUTC_TYPE) ,  allocatable,  dimension (:  )  ::  OUTPUTC 

INTEGER  : :  ndata 

REAL,  PARAMETER  ::  C_cutoff  =  l.E-20 
END  MODULE  datatype 
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in  which  (il,jl,kl)  and  (ig,jg,kg)  represent  (I,J,K)  indices  in  the  local  and  global 
coordinate  systems,  respectively.  Note  that  the  local  indices  (il,  jl,kl)  are  the  indices  on  a 
partition  (or  sub-domain)  of  the  solution  domain  designated  by  its  node  (or  CPU)  number:  node 
when  the  MPI  library  is  used  for  message  passing  on  a  parallel  Linux  cluster  in  SHARCNET. 
The  global  indices  (ig,jg,kg)  correspond  to  the  indices  on  a  single  grid  before  grid 
partitioning  is  done  by  the  pre-processor:  urbanPartitioning.  f  90 . 

A  header  file:  compact. h  listed  below  using  the  data  type  output and  outputc(:) 
defined  above  is  included  in  the  main  program  of  urbanSTREAM-P: 

ndata=0 

DO  K=NKBEG (NB) +  2 , NKEND (NB) -2 
DO  J=NJBEG (NB) +  2 , NJEND (NB) -2 
DO  I=NIBEG (NB) +  2 , NIEND (NB) -2 

C_dimensional=abs (C (I, J, K) ) *C_SOURCE  ! [kg/mA3] 

IF (C_dimensional<C_cutof f )  CYCLE 
ndata=ndata+l 
END  DO 
END  DO 
END  DO 

allocate (outputc (ndata) ) 

i 

i j  k_max= 

1  (NIEND (NB) -NIBEG (NB) -3) 

1* (NJEND (NB) -NJBEG (NB) -3) 

1* (NKEND (NB) -NKBEG (NB) -3) 

i 

n  =0 
nn=0 

DO  K=NKBEG (NB) +2 , NKEND (NB) -2 
DO  J=NJBEG (NB) +2 , NJEND (NB) -2 
DO  I=NIBEG (NB) +2 , NIEND (NB) -2 
nn=nn+l 

C_dimensional=abs (C(I, J,K) ) *C_SOURCE  ! [kg/mA3] 

IF (C_dimensional<C_cutof f )  CYCLE 
n=n+l 

outputc (n) %i=output (nn, node+1 ) %IG 
outputc (n) % j  =output (nn, node+1 )  % JG 
outputc (n) %k=output (nn, node+1 ) %KG 
outputc (n) %c=C_dimensional 
END  DO 
END  DO 
END  DO 

i 

write (30)  nint ( itime*delt*Tscal ) 
write  (30)  ndata 
if(ndata/=0)  then 

write (30)  (outputc (n) %i, outputc (n) %j , outputc (n) %k,  outputc (n) %c, 
ln=l , ndata) 
end  if 

i 

deallocate (outputc) 


compute  total  number  of 
concentration  data  where 
C>C_cutoff  is  exported 
to  urbanCollect . f 90  for 
post-processing 
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Note  that  the  info  about  global  indices  (ig,jg,kg)  are  provided  by  the  pre-processor: 
urbanPartitioning. f 90  through  the  subroutine  geometry. f 90  in  urbanSTREAM-P  shown 
below: 


read(io) n, nnmax 

allocate (output (nnmax, nbloc)  ) 

read ( io) ( ( 

loutput  (nn, n) %IL, output (nn, n) % JL, output (nn, n) %KL, 
loutput  (nn, n) %IG, output (nn, n) % JG, output (nn, n) %KG) , 
lnn=l , nnmax) 


The  output  of  concentration  fields  obtained  with  urbanSTREAM-P  from  each  CPU  (or  node  in 
the  sample  code)  is  collected  by  urbanCoiiect.f90  to  map  (or  merge)  solutions  from  each 
partition  onto  a  single  mesh  of  498x336x65  grid  points  before  post-processing  the  results  using 
Tecplot.  The  main  part  of  urbanCoiiect .  f  90  is  listed  below  for  illustration  purpose: 


C 


C 


C 


C 


read (30+node)  time 
read (30+node)  ndata 


if(ndata/=0)  then 
allocate (outputc (ndata) ) 


from  urbanSTREAM-P 

— o - 


read (30+node) (outputc (n) %i, outputc (n) %j , outputc (n) %k, outputc (n) %c, 
ln=l , ndata) 


do  n=l, ndata 
IG=outputc (n)  %i 
JG=outputc (n)  %j 
KG=outputc (n) %k 
C ( IG, JG, KG) =outputc (n) %c 
end  do 


deallocate (outputc) 
end  if 


c 

end  do  lend  of  block-loop 

i 


!  tecplot  format 

i 


i0= (time -mod (time,  1000)  )  /1000 

il= (time-iO* 1000-mod (time, 100 ) ) /100 

i2  = ( time-iO* 1000-il* 100-mod ( time-iO* 1000-il* 100, 10) ) /10 
i3=time-i0* 1000-il* 100-i2* 10 


!  comp 
!  comp 
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END  DO 
END  DO 
END  DO 
close  ( 60 ) 


In  summary,  the  compact  data  structure  developed  in  the  present  project  is  added  to 
urbanSTREAM-P,  the  parallel  version  of  urbanSTREAM,  urbanPartitioning.f90  and 
urbanCollect .  f  90 .  As  shown  in  Figure  79,  urbanPartitioning.f90  and  urban - 
Collect .  f  90  are  pre-processor  and  post-processor  of  urbanSTREAM-P  interfacing  urbanGRID 
and  urbanPOST. 


4.4.2  Validation  and  storage  efficiency  test 

Sample  iso-surfaces  of  the  concentration  field  at  log(C)  =  -3,  -5  and  -8  for  a  southwest  (or  225°) 
wind  direction  for  Scenario  1  are  shown  in  Figures  80  and  81  at  t  =  1,800  sec,  where  contours  of 
log(C)  below  -12  are  not  shown  as  they  are  insignificant  for  practical  applications.  As  seen  from 
Figure  80  (no  compact  data  structure  is  used)  and  Figure  81  (the  proposed  compact  data  structure 
is  adopted),  the  results  are  essentially  identical.  The  time  history  of  the  compression  ratio  is 
given  in  Figure  82.  When  time  is  small,  say,  t  =  60  sec,  the  compression  ratio  is  as  low  as  0.8%. 
As  the  time  approaches  to  t  =  1,800  sec,  the  compression  ratio  reaches  about  27%. 


4.4.3  File  transfer  and  dataset  naming  conventions 

In  order  to  transfer  (via  ftp)  the  dataset  for  each  scenario  at  different  times  up  to  t  =  40  min  to 
3DIntemet  Inc.,  a  Fortran  code:  read_data .  f  90  was  provided  to  3DIntemet  Inc.  to  reproduce 
our  results  generated  by  Tecplot.  The  key  portion  of  the  code  is  listed  below  to  convert  the 
concentration  data  from  the  compact  data  structure  to  the  standard  (i ,  j,k)  array  data  structure: 


do  n=l,ndata 

C (output (n) %i, output (n) %j , output (n) %k) ^output  (n) %c 
end  do 


The  complete  listing  of  the  Fortran  program  is  given  in  Appendix  A  for  information  only. 
Sample  concentration  results  obtained  from  urbanPOST  and  read_data.f90,  respectively,  are 
given  in  Figures  83,  and  identical  results  are  seen  to  ensure  the  validity  of  the  code. 

For  each  scenario  with  a  given  source  location,  two  sets  of  data  are  provided  individually; 
namely,  one  for  the  continuous  release  and  the  other  for  the  puff  release  as  explained  below.  For 
Scenario  1  discussed  in  Section  3.1,  the  release  can  be  split  into  two  parts: 

•  Release  1  (Rl):  30,000  kg  over  30  min  or  1,000  kg/min  (=16.667  kg/s)  over  30  min  at  At 
=  10  sec,  which  is  considered  as  “continuous  release”. 

•  Release  2  (R2):  30,000  kg  over  one  time  step:  At  =  10  sec,  which  is  considered  as 
“instantaneous  puff  release”  or  “puff  release”. 

Since  the  concentration  equation  is  a  linear  PDE  (partial-differential  equation),  with  a  given 
velocity  field  the  final  solution  of  the  concentration  field  is  obtained  by  combining  Rl  and  R2 
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releases  together.  Iso-surfaces  of  the  concentration  field  at  logio(C)=  -2.2,  -3.3  and  -4.4  for  a 
source  of  continuous  release,  of  puff  release  and  of  combining  continuous  and  puff  releases 
together  at  the  northwest  wind  direction  are  shown  in  Figures  84,  85  and  86,  respectively. 

Two  zip  files  are  produced  for  Scenario  1:  Sl-Rl-315.tar.gz  and  Sl-R2-315.tar.gz,  where  the 
naming  conventions  are: 

•  SI:  Scenario  1; 

•  Rl,  R2:  continuous  and  puff  releases,  respectively; 

•  315:315°  wind  direction. 

Similar  naming  conventions  apply  to  Scenario  2  at  6  other  wind  directions:  45°,  90°,  135°,  225°, 
270°  and  315°,  except  Rl,  R2  are  defined  differently  as  follows: 

•  Release  1  (Rl):  2  .7778  kg/s  over  30  min  at  At  =  10  sec,  which  is  considered  as 
“continuous  release”. 

•  Release  2  (R2):  10,000  kg  over  100  sec  (or  10  time  steps  at  At  =  10  sec),  which  is 
considered  as  “instantaneous  release”. 

For  each  release  (Rl  or  R2)  in  Scenarios  1  and  2  at  various  wind  directions,  270  data  files  in  total 
at  every  10  sec  for  0  <  t  <45  min  are  generated. 

For  Scenario  3,  concentration  results  with  only  one  wind  direction  at  135°  are  obtained  for 
source  locations  at  L3-1,  L3-2,  L3-3  and  L3-4  defined  in  Section  3.3.  Therefore,  datasets  named 
by  S3-Rl-L3-l.tar.gz  and  S3-Rl-L3-l.tar.gz  mean  that 

•  S3:  Scenario  3  at  southeast  (or  135°)  wind  direction; 

•  Rl,  R2:  continuous  and  puff  releases,  respectively; 

•  L3-1:  source  located  at  L3-1  [or  at  (51  deg,  02’,  39.18”  N,  114  deg,  03’,  17.20”  W)]. 

Similarly,  for  source  locations  at  L4-1,  L4-2,  L4-3,  L4-4  and  L4-LC  defined  in  Section  3.3,  only 
one  wind  direction  at  225°  are  considered  to  generate  the  concentration  field  results.  Datasets, 
such  as,  S3-Rl-L4-l.tar.gz  and  S3-R1-L4-1  mean  that 

•  S3:  Scenario  3  at  southwest  (or  225°)  wind  direction; 

•  Rl,  R2:  continuous  and  puff  releases,  respectively; 

•  L4-1:  source  located  at  L4-1  [or  at  (51  deg,  02’,  41.46”  N,  114  deg,  05’,  21.06”  W)  ]. 

Rl  and  R2  above  for  Scenario  3  are  defined  as  follows: 

•  Release  1  (Rl):  50  kg  over  30  min  or  1.667  kg/min  (=0.02778  kg/s)  at  At  =  1  sec,  which 
is  considered  as  “continuous  release”. 

•  Release  2  (R2):  300  kg  over  one  time  step:  At  =  1  sec,  which  is  considered  as  “puff 
release”. 
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In  total,  348  data  files  for  each  release  (R1  or  R2)  in  Scenario  3  at  various  source  locations  are 
generated  according  to  the  following  frequency: 

•  every  1  sec  for  0<  t  <  1  min; 

•  every  5  sec  for  1  min  <  t  <  10  min; 

•  every  10  sec  for  10  min  <  t  <  40  min. 


5.  CONCLUSIONS  AND  RECOMMENDATIONS 

In  this  contract,  extensive  hazardous  plume  entities  have  been  computed  by  coupling 
urbanSTREAM  with  urbanEU  in  order  to  provide  the  resulting  concentration  datasets  over  a 
period  of  about  45  min  to  the  game-based  simulation  engine  developed  by  3DIntemet  Inc.  as  part 
of  the  CRTI  project  09-509TD.  urbanSTREAM  was  executed  in  a  standalone  mode.  In 
consequence,  at  the  inflow  (inlet)  boundary  of  the  computational  domain  a  power-law  velocity 
profile  described  in  Section  2.2  was  used  to  define  the  inflow  boundary  conditions  for  the 
simulations.  The  new  dense-gas  model  presented  in  Section  2  and  results  discussed  in  Section  4 
allow  the  following  conclusions  to  be  drawn,  and  recommendation  for  the  future  work  to  be 
made: 

1.  The  development  of  a  realistic  multi-spectrum  CBRNE  scenario  for  an  actual  urban 
environment  (downtown  Calgary)  involves  three  sequential  events  consisting  of  accidental 
and  deliberate  releases  of  hazardous  materials  (toxic  industrial  chemicals  and  chemical 
agent)  and  detonation  of  an  explosive  device,  which  are  referred  to  as  Scenarios  1  to  3  in 

Section  3.  There  are  42  concentration  files  (each  =10  GB  in  size)  that  need  to  be 

12 

transferred  to  3DIntemet  Inc.  via  FileTransfer  Protocol  (ftp)  utility  program.  Since  C<10‘ 
kg/m3  is  insignificant  for  practical  applications,  a  “compact  data  structure”  described  in 
Section  4.4  is  developed  to  facilitate  the  transfer  of  the  datasets.  For  example,  the 
compression  ratio  of  the  dataset  over  a  period  of  30  min  is  less  than  27%  (see  Figure  82) 
when  the  compact  data  structure  is  included  in  urbanSTREAM-P  [parallel  version  of 
urbanSTREAM  (Yee  et  al.,  2007)]  and,  its  preprocessor  and  postprocessors:  urban- 
partitioning,  f  90  and  urbanCollect . f 90. 

2.  For  the  dense-gas  model,  a  simple  and  rigorous  derivation  for  the  modeling  of  the 
buoyancy  term  within  the  framework  of  the  Boussinesq  approximation  is  provided  and 
results  obtained  with  this  model  has  been  published  in  the  International  Journal  of 
Greenhouse  Gas  Control  (Hsieh  et  al.,  2013).  The  effects  of  buildings  and/or  complex 
topography  on  the  dispersion  of  dense  gas  clouds/plumes  are  studied.  It  is  shown  that  the 
model  is  able  to  provide  fair  to  good  predictions  for  the  Thomey  Island  test  gas,  where 
experimental  data  is  available. 

3.  In  order  to  verify  the  implementation  of  dense-gas  model  in  urbanEU,  two  methods, 
referred  to  as  the  “Q-approach”  and  “Wjet-appraoch”,  are  attempted  here.  Preliminary 
dispersion  results  obtained  with  both  dense-gas  and  passive  gas  models  for  Scenario  1  in 
the  City  of  Calgary  are  fairly  comparable.  With  a  pre-computed  quasi-steady  wind  and 
turbulence  fields  as  the  initial  condition,  the  dense-gas  model  is  approximately  400  times 
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more  costly  than  the  passive-gas  approach.  For  the  latter,  the  advection-diffusion  equation 
for  the  passive  scalar  (or  concentration)  is  entirely  decoupled  from  the  other  6  equations 
(namely,  pressure-correction  equations,  3  momentum  equations,  k-  and  s-  equations  for 
turbulence).  In  other  words,  only  the  linear  advection-diffusion  equation  is  required  to  be 
solved  for  the  passive-gas  approach.  Therefore,  most  dispersion  results  presented  in  this 
report  are  obtained  with  the  passive-gas  approach. 

4.  To  further  validate  the  proposed  dense-gas  model,  LES  (large  eddy  simulations)  or  hybrid 
RANS/LES,  such  as  DES  (Detached  Eddy  Simulation)  [see,  e.g.,  Strelets  (2001)],  PRNS 
(Partially  Resolved  Numerical  Simulation)  [see,  e.g.,  Shih  and  Liu  (2008),  Hsieh  et  al. 
(2010)]  might  be  a  better  approach  than  the  present  (U)RANS  method.  However, 
application  of  LES  or  hybrid  RANS/LES  to  simulate  the  highly  disturbed  flow  in  a 
cityscape,  in  which  a  large  portion  of  buildings  is  explicitly  resolved,  can  be  prohibitively 
expensive  and  therefore  was  not  considered  in  this  project. 

5.  In  addition  to  the  Thomey  Island  field  trials,  which  was  conducted  from  July  1982  to  July 
1983,  more  recent  dense-gas  field  experiment,  such  as  the  industrial  scale  pipeline  rupture 
experiments  in  the  European  project  entitled  “CCfPipeHaz”  (http://www.co2pipehaz.eu/) 
[see,  also  Woolley  et  al.  (2013)],  should  be  considered  for  future  investigation  when  the 
data  becomes  publicly  available. 

6.  In  Scenario  2,  it  was  assumed  that  derailment  caused  10,000  kg  of  chlorine  to  be  expelled 
into  the  atmosphere  and  “flashed”  in  the  form  of  vapor/aerosol  droplets.  It  is  recommended 
that  a  (more  rigorous)  chlorine  source  term  model  [see,  e.g.,  Barrett  (2009)]  be  included  in 
urbanEU. 

7.  Scenario  3  involved  the  explosive  release  of  an  unknown  chemical  agent.  It  is 
recommended  that  an  unconfined  vapor  cloud  explosion  source  term  model  [see,  e.g.,  the 
SEVEX  project  (Ronday,  et  al.,  1995)]  be  developed  based  on  WATCFD’s  recent 
capability  in  CFD  modeling  of  various  explosion  problems  [see,  e.g.,  Ji  et  al.  (2010),  Xu  et 
al.  (2013)]  as  a  potential  extension  of  the  present  project. 
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Appendix  A:  Fortran  Code  Listing  for  read_data.f90 

MODULE  data 
!  »TYPE 

TYPE  OUTPUT_TYPE 
INTEGER  ::  I,  J,  K 
REAL : :  C 

END  TYPE  OUTPUT_TYPE 
!»TYPE 

END  MODULE  data 
PROGRAM  read_data 

i 

USE  data 

TYPE (OUTPUT_TYPE ) ,  allocatable,  dimension (: )  ::  OUTPUT 

INTEGER  : :  ndata 

i 

real,  dimension (:,:,:) ,  ALLOCATABLE : :  C 

real,  dimension ,  ALLOCATABLE::  X,Y,Z 

character*80  str 

open  ( 10 , f ile= ’ GRID . dat ’ ) 

read (10, *) NI,NJ,NK 

i 

ALLOCATE (X (NI,NJ,NK) , Y (NI,NJ,NK) , Z (NI,NJ,NK) ) 

ALLOCATE (C (NI,NJ,NK) ) 

i 

read (10,*) ( ( (X(I, J,K) , 1=1, NI) , J=1,NJ) ,K=1,NK) 
read (10,*) ( ( (Y(I, J,K) , 1=1, NI) , J=1,NJ) ,K=1,NK) 
read (10,*) ( ( (Z (I, J,K) ,I*1,NI) , J=1,NJ) ,K=1,NK) 

i 

CALL  GETARG (1, str) 

WRITE (*,*)  'Input  file  name:', str 
WRITE (*, *) 

OPEN (UNIT  =  11,  FILE  =  TRIM (str) ,  STATUS  =  'OLD',  FORM= ' FORMATTED ' , 

ACTION  =  'READ') 

read ( 11 , * )  ndata 
allocate (output (ndata) ) 

read (11, 20) (output (n) %i, output (n) % j , output (n) %k, output (n) %c, n=l, ndata) 

20  f ormat ( 315 , E15 . 8 ) 

close  ( 11 ) 

i 

write  (*,*) 'number  of  nozero  C>C_cutoff,  number  of  total  data,  & 
compression  ratio' 

write ( * ,  * ) ndata,  (ni ) * (n j ) * (nk) , & 

float (ndata) /float (  (ni) * (n j ) * (nk)  ) 

C=0  . 

do  n=l, ndata 

C  (output  (n) %i, output (n) %j , output (n) %k) =output (n) %c 
end  do 

i 

!  tecplot  format 

i 

open ( 60 , f ile= ' comp_' //str) 

WRITE (60,*) 'VARIABLES  =  "X" , "Y" , "Z" , "C" ' 
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WRITE (60,*) ’ ZONE  F=POINT,  I=',NI,  J=',NJ,  ',  K= ’ , NK 

DO  K=1 , NK 
DO  J=1 ,NJ 
DO  1=1, NI 

WRITE  (60,  *)X(I,  J,K)  ,  Y(I,  J,K)  ,  Z  (I,  J,K)  ,C(I,  J,K) 

END  DO 
END  DO 
END  DO 

END  PROGRAM  read  data 
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Table  1.  Parameters  used  for  the  power- law  specification  of  the  vertical  profile  of  the 
streamwise  mean  wind  velocity  for  Thomey  Island  trial  numbers  26  and  29  (Phase  II). 


Trial  number 

Wind  speed  u0  (m/s) 

Exponent  n 

26 

1.9 

0.07 

29 

5.6 

0.15 

32 
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Figure  1.  A  two-dimensional  view  of  the  computational  mesh  for  trial  number  26  in  a  vertical  x- 
y  plane  at  z  =  0  (top  panel)  and  in  a  horizontal  x-z  plane  at  y  =  0  (bottom  panel).  The  mesh 
consists  of  189  x  60  x  83  nodes  in  the  streamwise  (x),  vertical  (y),  and  spanwise  (z)  directions, 
respectively. 
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Trial  26:  t  =  0  s 
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Trial  26:  t  =  14.68  s 
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Trial  26:  t  =  57.32  s 
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Figure  2.  Time  evolution  of  mean  flow  streamlines  and  contours  of  mean  concentration  (%  v/v) 
of  the  gas  cloud  in  a  vertical  plane  at  z!H=  0  for  trial  number  26. 
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Trial  29:  t  =  0  s 
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Trial  29:  t  =  3.09  s 
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Trial  29:  t  =  8.47  s 
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Figure  3.  Time  evolution  of  mean  flow  streamlines  and  contours  of  mean  concentration  (%  v/v) 
of  the  gas  cloud  in  a  vertical  plane  at  z!H=  0  for  trial  number  29. 
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Figure  4.  Time  history  of  gas  concentration  on  the  windward  face  of  the  obstacle  at  a  height  of 
6.4  m  (top  panel)  and  on  the  leeward  face  of  the  obstacle  at  a  height  of  0.4  m  (bottom  panel)  for 
trial  number  26:  (-<>-)  experimental  data  from  Davies  and  Singh  (1985);  ( — )  simulation. 


36 


t(s) 

Figure  5.  Time  history  of  gas  concentration  on  the  leeward  face  of  the  obstacle  at  a  height  of  0.4 
m  for  trial  number  29:  (-<>-)  experimental  data  from  Davies  and  Singh  (1985);  ( — )  simulation. 
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Figure  6.  Buoyancy  production  term  sensitivity  analysis:  contours  of  turbulence  kinetic  energy 
(top  two  panels)  and  mean  concentration  (%  v/v)  in  the  gas  cloud  (bottom  two  panels)  in  a 
vertical  plane  at  z/H  =  0  for  trial  number  29,  with  and  without  inclusion  of  the  Gk  term  in  the 
turbulence  transport  equations. 
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Figure  7.  Google  map  showing  “inner  region”  and  “outer  region”  of  the  computational  domain. 


inner  region 


Figure  8.  Top  view  of  the  computational  grid  consisting  of  “inner  region”  and  “outer  region”. 
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Figure  9.  New  Shapefiles  including  the  Bow  Building  in  the  City  of  Calgary  provided  by 
3DIntemet  Inc. 


Figure  10.  3D  view  of  the  computational  grid  consisting  of  498x336x65  nodes  in  the  x-,  y-  and 
z-directions,  respectively. 
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Figure  11.  Sample  streamline  patterns  and  associated  velocity  fields  at  135°  wind  direction  in 
the  City  of  Calgary. 


Figure  12.  Sample  streamline  patterns  and  associated  velocity  fields  at  225°  wind  direction  in 
the  City  of  Calgary. 
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Figure  13.  Scenario  1  (front  view)  -  Sample  streamline  patterns  superimposed  with  iso-surfaces 
of  concentration  at  225°  wind  direction  in  the  City  of  Calgary. 


Figure  14.  Scenario  1  (rear  view)  -  Sample  streamline  patterns  superimposed  with  iso-surfaces 
of  concentration  at  225°  wind  direction  in  the  City  of  Calgary. 
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Figure  15.  Scenario  1  -  Contours  of  concentration  on  a  logarithmic  scale  at  t  =  5  min  at  a  wind 
direction  of  45°  with  a  source  located  at  (51  deg,  02’,  37.12”  N,  1 14  deg,  03’,  1 1.89”  W). 


Figure  16.  Scenario  1  -  Contours  of  concentration  on  a  logarithmic  scale  at  t  =  30  min  at  a  wind 
direction  of  45°  with  a  source  located  at  (51  deg,  02’,  37.12”  N,  1 14  deg,  03’,  1 1.89”  W). 
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Figure  17.  Scenario  1  -  Contours  of  concentration  on  a  logarithmic  scale  at  t  =  5  min  at  a  wind 
direction  of  90°  with  a  source  located  at  (51  deg,  02’,  37.12”  N,  1 14  deg,  03’,  1 1.89”  W). 


Figure  18.  Scenario  1  -  Contours  of  concentration  on  a  logarithmic  scale  at  t  =  30  min  at  a  wind 
direction  of  90°  with  a  source  located  at  (51  deg,  02’,  37.12”  N,  114  deg,  03’,  11.89”  W). 
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Figure  19.  Scenario  1  -  Contours  of  concentration  on  a  logarithmic  scale  at  t  =  5  min  at  a  wind 
direction  of  135°  with  a  source  located  at  (51  deg,  02’,  37.12”  N,  1 14  deg,  03’,  1 1.89”  W). 


Figure  20.  Scenario  1  -  Contours  of  concentration  on  a  logarithmic  scale  at  t  =  30  min  at  a  wind 
direction  of  135°  with  a  source  located  at(51  deg,  02’, 37. 12”  N,  114deg,  03’,  11.89”  W). 
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Figure  21.  Scenario  1  -  Contours  of  concentration  on  a  logarithmic  scale  at  t  =  5  min  at  a  wind 
direction  of  225°  with  a  source  located  at  (51  deg,  02’,  41.07”  N,  1 14  deg,  05’,  14.09”  W). 


Figure  22.  Scenario  1  -  Contours  of  concentration  on  a  logarithmic  scale  at  t  =  30  min  at  a  wind 
direction  of  225°  with  a  source  located  at  (51  deg,  02’,  41.07”  N,  114  deg,  05’,  14.09”  W). 
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Figure  23.  Scenario  1  -  Contours  of  concentration  on  a  logarithmic  scale  at  t  =  5  min  at  a  wind 
direction  of  270°  with  a  source  located  at  (51  deg,  02’,  41.07”  N,  1 14  deg,  05’,  14.09”  W). 


Figure  24.  Scenario  1  -  Contours  of  concentration  on  a  logarithmic  scale  at  t  =  30  min  at  a  wind 
direction  of270°  with  a  source  located  at(51  deg,  02’,  41.07”  N,  114deg,  05’,  14.09”  W). 
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Figure  25.  Scenario  1  -  Contours  of  concentration  on  a  logarithmic  scale  at  t  =  5  min  at  a  wind 
direction  of  315°  with  a  source  located  at  (51  deg,  02’,  41.07”  N,  1 14  deg,  05’,  14.09”  W). 


Figure  26.  Scenario  1  -  Contours  of  concentration  on  a  logarithmic  scale  at  t  =  30  min  at  a  wind 
direction  of  315°  with  a  source  located  at  (51  deg,  02’,  41.07”  N,  114  deg,  05’,  14.09”  W). 
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Figure  27.  Scenario  1  -  Contours  of  concentration  of  a  dense  gas  (chlorine)  at  t  =  10  min  at  a 
wind  direction  of  135°  with  a  source  at  Location  1  using  the  “Q-approach”. 


Figure  28.  Scenario  1  -  Contours  of  concentration  of  a  dense  gas  (chlorine)  at  t  =  10  min  at  a 
wind  direction  of  135°  with  a  source  at  Location  1  using  the  “Wjet-approach”. 
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Figure  29.  Scenario  1  (dense  gas  result)  -  Contours  of  log(C)  superimposed  with  the  velocity 
vector  field  at  t  =  10  min  at  a  wind  direction  of  135°  with  a  source  at  Location  1. 


Figure  30.  Scenario  1  (passive  gas  result)  -  Contours  of  log(C)  superimposed  with  the  velocity 
vector  field  at  t  =  10  min  at  a  wind  direction  of  135°  with  a  source  at  Location  1. 
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Figure  31.  Scenario  2  -  Contours  of  concentration  on  a  logarithmic  scale  at  t  =  5  min  at  a  wind 
direction  of  45°  with  a  source  located  at  (51  deg,  02’,  37.12”  N,  1 14  deg,  03’,  1 1.89”  W). 


Figure  32.  Scenario  2  -  Contours  of  concentration  on  a  logarithmic  scale  at  t  =  30  min  at  a  wind 
direction  of  45°  with  a  source  located  at  (51  deg,  02’,  37.12”  N,  1 14  deg,  03’,  1 1.89”  W). 
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Figure  33.  Scenario  2  -  Contours  of  concentration  on  a  logarithmic  scale  at  t  =  5  min  at  a  wind 
direction  of  90°  with  a  source  located  at  (51  deg,  02’,  37.12”  N,  1 14  deg,  03’,  1 1.89”  W). 


Figure  34.  Scenario  2  -  Contours  of  concentration  on  a  logarithmic  scale  at  t  =  30  min  at  a  wind 
direction  of  90°  with  a  source  located  at  (51  deg,  02’,  37.12”  N,  114  deg,  03’,  11.89”  W). 
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Figure  35.  Scenario  2  -  Contours  of  concentration  on  a  logarithmic  scale  at  t  =  5  min  at  a  wind 
direction  of  135°  with  a  source  located  at  (51  deg,  02’,  37.12”  N,  1 14  deg,  03’,  1 1.89”  W). 


Figure  36.  Scenario  2  -  Contours  of  concentration  on  a  logarithmic  scale  at  t  =  30  min  at  a  wind 
direction  of  135°  with  a  source  located  at  (51  deg,  02’,  37.12”  N,  1 14  deg,  03’,  1 1.89”  W). 
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Figure  37.  Scenario  2  -  Contours  of  concentration  on  a  logarithmic  scale  at  t  =  5  min  at  a  wind 
direction  of  225°  with  a  source  located  at  (51  deg,  02’,  41.07”  N,  1 14  deg,  05’,  14.09”  W). 


Figure  38.  Scenario  2  -  Contours  of  concentration  on  a  logarithmic  scale  at  t  =  30  min  at  a  wind 
direction  of  225°  with  a  source  located  at  (51  deg,  02’,  41.07”  N,  114  deg,  05’,  14.09”  W). 
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Figure  39.  Scenario  2  -  Contours  of  concentration  on  a  logarithmic  scale  at  t  =  5  min  at  a  wind 
direction  of  270°  with  a  source  located  at  (51  deg,  02’,  41.07”  N,  1 14  deg,  05’,  14.09”  W). 


Figure  40.  Scenario  2  -  Contours  of  concentration  on  a  logarithmic  scale  at  t  =  30  min  at  a  wind 
direction  of  270°  with  a  source  located  at  (51  deg,  02’,  41.07”  N,  114  deg,  05’,  14.09”  W). 
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Figure  41.  Scenario  2  -  Contours  of  concentration  on  a  logarithmic  scale  at  t  =  5  min  at  a  wind 
direction  of  315°  with  a  source  located  at  (51  deg,  02’,  41.07”  N,  1 14  deg,  05’,  14.09”  W). 


Figure  42.  Scenario  2  -  Contours  of  concentration  on  a  logarithmic  scale  at  t  =  30  min  at  a  wind 
direction  of  315°  with  a  source  located  at  (51  deg,  02’,  41.07”  N,  114  deg,  05’,  14.09”  W). 
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Datasets  delivered  to  3Dlnternet  [At=l  sec] 


51 0442166667  N  1 14  0547777778  W 


Wind  direg 
at  L3-lr 


*n:  135  degrees 


^  III:  Van  containing  two  50-US  gallon  drums  of  ^ 
a  chemical  warfare  agent  or  CWA  (e.g.,  sarin, 
mustard  etc.) 

\ _ _ _ / 


Figure  43.  Scenario  3  -  Location  of  the  source  L3-1  at  (51  deg,  02’,  39.18”  N,  114  deg,  03’, 
17.20”  W)  at  a  wind  direction  of  135°. 


Figure  44.  Scenario  3  -  Top  view  of  contours  of  concentration  on  a  logarithmic  scale  at  t  =  5 
min  at  a  wind  direction  of  135°  with  the  source  L3-1  located  at  (51  deg,  02’,  39.18”  N,  1 14  deg, 
03’,  17.20”  W). 
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Figure  45.  Scenario  3  -  Contours  of  concentration  on  a  logarithmic  scale  at  t  =  5  min  at  a  wind 
direction  of  135°  with  the  source  L3-1  located  at  (51  deg,  02’,  39.18”  N,  114  deg,  03’,  17.20” 
W). 


Figure  46.  Scenario  3  -  Contours  of  concentration  on  a  logarithmic  scale  at  t  =  15  min  at  a  wind 
direction  of  135°  with  the  source  L3-1  located  at  (51  deg,  02’,  39.18”  N,  114  deg,  03’,  17.20” 
W). 
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Datasets  delivered  to  3Dlnternet  [At=l  sec] 


51 0440611111  N  114  0556056656  W 


^  III:  Van  containing  two  50-US  gallon  drums  of  ^ 
a  chemical  warfare  agent  or  CWA  (e.g.,  sarin, 
mustard  etc.) 

\ _ _ _ / 


Figure  47.  Scenario  3  -  Location  of  the  source  L3-2  at  (51  deg,  02’,  38.62”  N,  114  deg,  03’, 
20.18”  W)  at  a  wind  direction  of  135°. 
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Figure  48.  Scenario  3  -  Top  view  of  contours  of  concentration  on  a  logarithmic  scale  at  t  =  5 
min  at  a  wind  direction  of  135°  with  the  source  L3-2  located  at  (51  deg,  02’,  38.62”  N,  1 14  deg, 
03’,  20. 18”  W). 
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Figure  49.  Scenario  3  -  Contours  of  concentration  on  a  logarithmic  scale  at  t  =  5  min  at  a  wind 
direction  of  135°  with  the  source  L3-2  located  at  (51  deg,  02’,  38.62”  N,  114  deg,  03’,  20.18” 
W). 


Figure  50.  Scenario  3  -  Contours  of  concentration  on  a  logarithmic  scale  at  t  =  15  min  at  a  wind 
direction  of  135°  with  the  source  L3-2  located  at  (51  deg,  02’,  38.62”  N,  114  deg,  03’,  20.18” 
W). 
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Datasets  delivered  to  3Dlnternet  [At=l  sec] 


Wind  direction:  135  degrees 
at  L3-3 


^  III:  Van  containing  two  50-US  gallon  drums  of  ^ 
a  chemical  warfare  agent  or  CWA  (e.g.,  sarin, 
mustard  etc.) 

\ _ _ _ / 


Figure  51.  Scenario  3  -  Location  of  the  source  L3-3  at  (51  deg,  02’,  35.87”  N,  114  deg,  03’, 
16.91”  W)  at  a  wind  direction  of  135°. 


Figure  52.  Scenario  3  -  Top  view  of  contours  of  concentration  on  a  logarithmic  scale  at  t  =  5 
min  at  a  wind  direction  of  135°  with  the  source  L3-3  located  at  (51  deg,  02’,  35.87”  N,  1 14  deg, 
03’,  16.91”  W). 
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Figure  53.  Scenario  3  -  Contours  of  concentration  on  a  logarithmic  scale  at  t  =  5  min  at  a  wind 
direction  of  135°  with  the  source  L3-3  located  at  (51  deg,  02’,  35.87”  N,  114  deg,  03’,  16.91” 
W). 


Figure  54.  Scenario  3  -  Contours  of  concentration  on  a  logarithmic  scale  at  t  =  15  min  at  a  wind 
direction  of  135°  with  the  source  L3-3  located  at  (51  deg,  02’,  35.87”  N,  114  deg,  03’,  16.91” 
W). 


62 


Datasets  delivered  to  3Dlnternet  [At=l  sec] 


51 0432805556  N  1 14  0654333333  W 


Wind  directionr"135  degrees 
at  LB-4 


'  III:  Van  containing  two  50-US  gallon  drums  of  ' 
a  chemical  warfare  agent  or  CWA  (e.g.,  sarin, 
mustard  etc.) 

_ _ _ / 


Figure  55.  Scenario  3  -  Location  of  the  source  L3-4  at  (51  deg,  02’,  35.81”  N,  114  deg,  03’, 
19.56”  W)  at  a  wind  direction  of  135°. 


Figure  56.  Scenario  3  -  Top  view  of  contours  of  concentration  on  a  logarithmic  scale  at  t  =  5 
min  at  a  wind  direction  of  135°  with  the  source  L3-4  located  at  (51  deg,  02’,  35.81”  N,  114  deg, 
03’,  19.56”  W). 
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Figure  57.  Scenario  3  -  Contours  of  concentration  on  a  logarithmic  scale  at  t  =  5  min  at  a  wind 
direction  of  225°  with  the  source  L3-4  located  at  (51  deg,  02’,  35.81”  N,  114  deg,  03’,  19.56” 
W). 


Figure  58.  Scenario  3  -  Contours  of  concentration  on  a  logarithmic  scale  at  t=  15  min  at  a  wind 
direction  of  225°  with  the  source  L3-4  located  at  (51  deg,  02’,  35.81”  N,  114  deg,  03’,  19.56” 
W). 
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Datasets  delivered  to  3Dlnternet  [At=l  sec] 


51  0448500000  H1 14  08018J3333  W 


Wind  direction:  225  degrees 
at  L4-1 

V _ 


III:  Van  containing  two  50-US  gallondrumsof^ 
a  chemical  warfare  agent  or  CWA  (e.g.,  sarin, 
mustard  etc.) 

\ _ _ _ / 


Figure  59.  Scenario  3  -  Location  of  the  source  L4-1  at  (51  deg,  02’,  41.46”  N,  114  deg,  05’, 
21.06”  W)  at  a  wind  direction  of  225°. 


Figure  60.  Scenario  3  -  Top  view  of  contours  of  concentration  on  a  logarithmic  scale  at  t  =  5 
min  at  a  wind  direction  of  225°  with  the  source  L4-1  located  at  (51  deg,  02’,  41.46”  N,  1 14  deg, 
05’,  21.06”  W). 
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Figure  61.  Scenario  3  -  Contours  of  concentration  on  a  logarithmic  scale  at  t  =  5  min  at  a  wind 
direction  of  225°  with  the  source  L4-1  located  at  (51  deg,  02’,  41.46”  N,  114  deg,  05’,  21.06” 
W). 


Figure  62.  Scenario  3  -  Contours  of  concentration  on  a  logarithmic  scale  at  t  =  15  min  at  a  wind 
direction  of  225°  with  the  source  L4-1  located  at  (51  deg,  02’,  41.46”  N,  114  deg,  05’,  21.06” 
W). 
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Datasets  delivered  to  3Dlnternet  [At=l  sec] 


51  044*055556  N  1 14  0882527778  W 


Wind  direjjtkSn:  225  degrees 
at  L4-2 

V _ 


III:  Van  containing  two  50-US  gallondrumsof^ 
a  chemical  warfare  agent  or  CWA  (e.g.,  sarin, 
mustard  etc.) 

\ _ _ _ ✓ 


Figure  63.  Scenario  3  -  Location  of  the  source  L4-2  at  (51  deg,  02’,  41.30”  N,  114  deg,  05’, 
17.71”  W)  at  a  wind  direction  of  225°. 


Figure  64.  Scenario  3  -  Top  view  of  contours  of  concentration  on  a  logarithmic  scale  at  t  =  5 
min  at  a  wind  direction  of  225°  with  the  source  L4-2  located  at  (51  deg,  02’,  41.30”  N,  1 14  deg, 
05’,  17.71”  W). 
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Figure  65.  Scenario  3  -  Contours  of  concentration  on  a  logarithmic  scale  at  t  =  5  min  at  a  wind 
direction  of  225°  with  the  source  L4-2  located  at  (51  deg,  02’,  41.30”  N,  114  deg,  05’,  17.71” 
W). 


Figure  66.  Scenario  3  -  Contours  of  concentration  on  a  logarithmic  scale  at  t  =  15  min  at  a  wind 
direction  of  225°  with  the  source  L4-2  located  at  (51  deg,  02’,  41.30”  N,  114  deg,  05’,  17.71” 
W). 
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Datasets  delivered  to  3Dlnternet  [At=l  sec] 


51 0446027778  N 


225  degrees 


^  III:  Van  containing  two  50-US  gallon  drums  of  ^ 
a  chemical  warfare  agent  or  CWA  (e.g.,  sarin, 
mustard  etc.) 

\ _ _ _ ✓ 


Figure  67.  Scenario  3  -  Location  of  the  source  L4-3  at  (51  deg,  02’,  40.57”  N,  114  deg,  05’, 
20.43”  W)  at  a  wind  direction  of  225°. 


Figure  68.  Scenario  3  -  Top  view  of  contours  of  concentration  on  a  logarithmic  scale  at  t  =  5 
min  at  a  wind  direction  of  225°  with  the  source  L4-3  located  at  (51  deg,  02’,  40.57”  N,  1 14  deg, 
05’,  20.43”  W). 
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Figure  69.  Scenario  3  -  Contours  of  concentration  on  a  logarithmic  scale  at  t  =  5  min  at  a  wind 
direction  of  225°  with  the  source  L4-3  located  at  (51  deg,  02’,  40.57”  N,  114  deg,  05’,  20.43” 
W). 


Figure  70.  Scenario  3  -  Contours  of  concentration  on  a  logarithmic  scale  at  t=  15  min  at  a  wind 
direction  of  225°  with  the  source  L4-3  located  at  (51  deg,  02’,  40.57”  N,  114  deg,  05’,  20.43” 
W). 
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Datasets  delivered  to  3Dlnternet  [At=l  sec] 


Wind  direction:  225  degrees 
at  L4-4 


^  III:  Van  containing  two  50-US  gallon  drums  of  ^ 
a  chemical  warfare  agent  or  CWA  (e.g.,  sarin, 
mustard  etc.) 

\ _ _ _ ✓ 


Figure  71.  Scenario  3  -  Location  of  the  source  L4-4  at  (51  deg,  02’,  40.50”  N,  114  deg,  05’, 
18.24”  W)  at  a  wind  direction  of  225°. 


Figure  72.  Scenario  3  -  Top  view  of  contours  of  concentration  on  a  logarithmic  scale  at  t  =  5 
min  at  a  wind  direction  of  225°  with  the  source  L4-4  located  at  (51  deg,  02’,  40.50”  N,  1 14  deg, 
05’,  18.24”  W). 
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Figure  73.  Scenario  3  -  Contours  of  concentration  on  a  logarithmic  scale  at  t  =  5  min  at  a  wind 
direction  of  225°  with  the  source  L4-4  located  at  (51  deg,  02’,  40.50”  N,  114  deg,  05’,  18.24” 
W). 


Figure  74.  Scenario  3  -  Contours  of  concentration  on  a  logarithmic  scale  at  t  =  15  min  at  a  wind 
direction  of  225°  with  the  source  L4-4  located  at  (51  deg,  02’,  40.50”  N,  114  deg,  05’,  18.24” 
W). 
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Datasets  delivered  to  3Dlnternet  [At=l  sec] 


Figure  75.  Scenario  3  -  Location  of  the  source  L4-LC  at  (5 1  deg,  02’,  40.96”  N,  1 14  deg,  05’, 
19.48”  W)  at  a  wind  direction  of  225°. 


Figure  76.  Scenario  3  -  Top  view  of  contours  of  concentration  on  a  logarithmic  scale  at  t  =  5 
min  at  a  wind  direction  of  225°  with  the  source  L4-LC  located  at  (51  deg,  02’,  40.96”  N,  114 
deg,  05’,  19.48”  W). 
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Figure  77.  Scenario  3  -  Contours  of  concentration  on  a  logarithmic  scale  at  t  =  5  min  at  a  wind 
direction  of  225°  with  the  source  L4-LC  located  at  (51  deg,  02’,  40.96”  N,  1 14  deg,  05’,  19.48” 
W). 


Figure  78.  Contours  of  concentration  on  a  logarithmic  scale  at  t  =  15  min  at  a  wind  direction  of 
225°  with  the  source  L4-LC  located  at  (51  deg,  02’,  40.96”  N,  1 14  deg,  05’,  19.48”  W). 
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!comp 

MODULE  datatype 
integer,  dimension^), 
ALLOCATABLE: :  NNM  AX 
!»TYPE 

TYPE  OUTPUT_TYPE 
INTEGER  ::  IL,  JL,  KL 
INTEGER  ::  IG,  JG,  KG 
END  TYPE  OUTPUT_TYPE 
!»TYPE 

END  MODULE  datatype 
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Figure  79.  Illustration  of  compact  data  structure  being  used  in  urbanPartitioning.f90, 
urbanCollect .  f  90  and  urbanSTREAM-P.  Note  that  urbanPartitioning.f90  and 
urbanCollect .  f 90  are  pre-processor  and  post-processor  of  urbanSTREAM-P  interfacing 
urbanGRID  and  urbanPOST. 


Figure  80.  Scenario  1  -  Contours  of  concentration  on  a  logarithmic  scale  at  t  =  30  min  for 
southwest  wind  direction.  No  compact  data  structure  is  used  (viz.,  compression  ratio  is  100%). 
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Figure  81.  Scenario  1  -  Contours  of  concentration  on  a  logarithmic  scale  at  t  =  30  min  for 
southwest  wind  direction.  The  compact  data  structure  is  employed,  and  the  compression  ratio  is 
27%  with  Ccut.off=l.E-20  kg/m3. 


Figure  82.  Time  history  of  compression  ratio  for  the  concentration  data  using  the  compact  data 
structure  for  Scenario  1  (see  also  Figure  81). 
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urbanSTREAM -P  ->  urbanCollect.f90  ->  urbanPOST 

1  I 


read  data.f90 


Identical 


Figure  83.  Scenario  1  -  Sample  concentration  results  obtained  from  urbanPOST  and 
read_data .  f  90  for  west  wind  direction  (or  270°  wind  direction)  at  t  =  300  sec  to  show  the 
validity  of  the  data  converter:  read_data .  f  90  (see  also  Appendix  A). 


Figure  84.  Scenario  1  -  Contours  of  concentration  on  a  logarithmic  scale  at  t  =  5  min  at  a  wind 
direction  of  315°  from  a  source  of  continuous  release. 
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Figure  85.  Scenario  1  -  Contours  of  concentration  on  a  logarithmic  scale  at  t  =  5  min  at  a  wind 
direction  of  3 15°  from  a  source  of  puff  release. 


Figure  86.  Scenario  1  -  Contours  of  concentration  on  a  logarithmic  scale  at  t  =  5  min  at  a  wind 
direction  of  3 15°  from  a  source  of  combining  both  continuous  and  puff  releases. 
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